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7.1  First  order  reconstructions  of  four  objects  at  the  limit  of  the 

Born  series  are  shown  here.  Objects  with  a  reconstruction  worse 
then  these  can  not  be  improved  by  an  inversion  procedure  based  on 
the  Born  series  because  the  Born  series  will  not  converge . 192 


ABSTRACT 


This  work  reviews  the  theory  and  limits  of  first  order  diffraction  tomogra¬ 
phy  and  studies  iterative  techniques  that  can  be  used  to  improve  the  quality  of 
tomographic  imaging  with  diffracting  sources.  Conventional  (straight-ray) 
tomographic  algorithms  are  not  valid  when  used  with  acoustic  or  microwave 
energy.  Thus  more  sophisticated  algorithms  are  needed. 

First  order  diffraction  tomography  uses  a  linearized  version  of  the  wave 
equation  and  gives  an  especially  simple  reconstruction  algorithm.  This  work 
reviews  first  order  approximations  to  the  scattered  field  and  studies  the  quality 
of  the  reconstructions  when  the  assumptions  behind  these  approximations  are 
violated.  It  will  be  shown  that  the  Bom  approximation  is  valid  when  the  phase 
change  across  the  object  is  less  than  t  and  the  Rytov  approximation  is  valid 
when  the  refractive  index  changes  by  less  than  two  or  three  percent. 

Better  reconstructions  will  be  based  on  higher  order  approximations  to  the 
scattered  field.  This  work  describes  two  fixed  point  algorithms  (the  Born  and 
the  Rytov  approximations)  and  an  algebraic  approach  to  more  accurately  cal¬ 
culate  the  scattered  fields.  The  limits  of  each  of  these  approaches  is  discussed 
and  simulated  results  are  shown. 

Finally  a  review  of  higher  order  inversion  techniques  is  presented.  Each  of 
these  techniques  is  reviewed  and  some  of  their  limitations  are  discussed. 


The  word  tomography  comes  from  the  Greek  words  tomo,  meaning 
sectional,  and  graphy,  meaning  representation.  Thus  a  tomographic  image  is  a 
cross  sectional  image  of  an  object.  As  the  term  is  used  today  tomography 
refers  to  a  procedure  to  collect  data  about  the  internal  structure  of  an  object 
and  then  mathematically  generate  an  image  of  some  otherwise  hidden  property 
of  the  object. 

Diffraction,  on  the  other  hand,  describes  the  spreading  of  acoustic  and 
electromagnetic  waves  as  they  propagate  through  space  and  around  objects. 
While  conventional  tomography  uses  x-rays  to  generate  an  image  of  the 
object’s  x-ray  attenuation  other  sources  of  energy  can  also  be  used.  Thus 
diffraction  tomography  uses  diffracting  energy  sources  to  illuminate  the  object 
and  then  generates  a  cross  sectional  image  of  the  object.  Since  ultrasound  and 
microwaves  diffract  and  refract  as  they  pass  through  most  objects  they  require 
more  sophisticated  algorithms  then  the  ones  used  for  x-ray  tomography.  These 
new  algorithms  for  diffraction  tomography  are  the  subject  of  this  work. 

Tomography  first  became  practical  only  a  few  years  ago  with  the  invention 
of  the  CAT  (Computer  Assisted  Tomography)  scanner  [Hou72|.  Hounsfield 
implemented  a  machine  that  illuminated  an  object  with  x-rays  and  measured 
the  proportion  of  energy  that  passed  through  the  object.  Then  by  inverting  a 
large  system  of  equations  he  was  able  to  generate  an  accurate  estimate  of  the 
spatial  variation  of  x-ray  attenuation  in  the  object. 

The  ability  to  generate  a  tomographic  image  of  an  object  has 
revolutionized  the  medical  field.  For  the  first  time  it  was  possible  to  get  a  clear 
image  of  the  internal  morphology  of  a  patient  without  the  use  of  surgery. 
Now,  x-ray  CAT  scanners  are  routinely  built  with  resolutions  of  less  than  a 
millimeter  and  images  with  more  than  512x512  pixels  [Kak85,  Her80,  Mac83, 
Bar81]. 

While  medical  CAT  scanners  often  generate  an  image  of  an  object’s  x-ray 
attenuation  there  are  limitations  to  this  procedure.  Foremost  is  the  fact  that 


not  all  types  of  soft  tissue  are  differentiated  by  their  x-ray  attenuation.  Thus 
x-ray  CAT  scans  have  wide  use  for  orthopedic  medicine  but  are  of  limited  use, 
for  example,  in  diagnosing  malignant  vs.  benign  tumors.  In  addition  x-rays  are 
an  ionizing  radiation  and  thus  there  is  a  small  chance  of  cancer  with  each  use. 
This  prevents,  for  example,  the  use  of  x-ray  CAT  scanners  for  mass  screening 
of  female  patients  for  breast  cancer. 

X-ray  tomography  is  based  on  the  Fourier  Slice  Theorem.  Consider  the 
experiment  shown  on  the  left  side  of  Figure  1.1.  Here  a  projection  is  shown 
that  represents  the  attenuation  of  the  object  along  each  of  the  indicated  lines. 
The  Fourier  Slice  Theorem  states  that  the  Fourier  transform  of  the  projection 
is  equal  to  the  values  of  the  two  dimensional  Fourier  transform  of  the  object 
along  a  radial  line.  An  estimate  of  the  object  can  then  be  formed  by  measuring 
projections  at  a  number  of  angles  and  then  simply  inverting  the  Fourier  data. 

Conventional  tomography  is  based  on  the  idea  that  x-rays  travel  in 
straight  lines  through  the  object  and  a  projection  measures  the  total  x-ray 
attenuation  of  the  object  along  straight  lines.  When  the  object  is  relatively 
large  and  has  a  small  refractive  index  it  is  possible  to  use  other  types  of  energy, 
for  example  microwaves,  seismic  and  ultrasound,  to  image  the  object.  With  a 
small  refractive  inex  the  energy  doesn’t  bend  as  it  goes  through  the  object  and 
thus  it  is  possible  to  measure  the  attenuation  of  the  object  along  straight  lines. 
This  is  the  only  requirement  needed  to  use  the  Fourier  Slice  Theorem  and  form 
an  image  of  the  object’s  acoustic  or  microwave  attenuation  [Gre74,  Gre75, 
Car76,  Jak76,  Glo77  and  Cra82j. 

Since  microwaves  and  acoustic  waves  are  easier  to  generate  and  measure 
than  x-rays  it  is  also  possible  to  generate  images  of  the  object’s  refractive 
index.  As  was  mentioned  earlier  it  is  necessary  to  assume  that  the  refractive 
index  change  is  small  so  that  the  energy  doesn’t  bend  as  it  travels  through  the 
object.  If  a  projection  is  measured  representing  the  delay  encountered  by  the 
energy  as  it  travels  through  different  parts  of  the  object  then  an  image  is 
formed  of  the  object’s  acoustic  or  electromagnetic  refractive  index.  This  extra 
information  can  often  make  it  easier  to  characterize  the  object. 

Two  methods  have  been  used  to  form  images  when  the  energy  no  longer 
travels  through  the  object  in  a  straight  line.  Perhaps  the  most  straightforward 
approach  is  to  simply  model  the  flow  of  energy  through  the  object  as  a  ray  and 
calculate  its  location  based  on  the  refractive  index  of  the  object  [And82,  Her76, 
Her73],  Unfortunately  these  algorithms  can  only  be  used  when  the  refractive 
index  change  is  less  than  10  or  20  percent  and  most  of  the  energy  is  refracted 
instead  of  diffracted.  Thus  this  approach  is  only  valid  when  the  wavelength  of 
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the  energy  is  much  smaller  than  any  details  of  the  object  [And84]. 

A  second  approach  is  to  model  the  flow  of  energy  through  the  object  with 
the  wave  equation.  While  this  approach  is  more  accurate  then  other  approaches 
it  is  not  always  possible  to  invert  the  resulting  system  of  equations  and  find  a 
closed  form  solution.  This  is  the  core  of  the  problem  for  successful  diffraction 
tomographic  images. 

A  simple  approach  to  solve  the  wave  equation  is  to  linearize  it  [Ish78, 
Che60,  Sla84,  Mue79,  Wol69j.  This  is  usually  done  by  assuming  that  the 
object  represents  a  very  small  perturbation  to  the  field.  Only  the  linear  terms 
are  retained  and  all  higher  order  terms  are  simply  ignored.  Unfortunately  this 
approach  is  also  limited  to  those  object  that  satisfy  the  constraints  of  the 
approximations.  As  will  be  shown  later  in  this  work  linearizing  the  wave 
equation  greatly  limits  the  objects  that  can  be  imaged. 

Finally  in  the  past  years  work  has  been  done  on  iterative  techniques  to 
solve  the  wave  equation.  Most  of  the  work  was  originally  applied  to  the 
inverse  scattering  problem  of  high  energy  physics  [Bal78,  New66,  Tay83]  and 
only  recently  applied  to  the  diffraction  tomography  problem. 

This  work  is  in  three  parts:  the  derivation  of  the  wave  equation  and  first 
order  reconstruction  algorithms,  limitations  of  first  order  algorithms  and  finally 
a  summary  of  iterative  techniques  that  can  be  applied  to  the  diffraction 
tomography  problem. 

First  in  Chapter  2  the  wave  equation  is  defined  for  both  acoustic  and 
electromagnetic  experiments.  This  scalar  equation  is  valid  for  both  types  of 
energy  and  forms  the  basis  of  all  work  to  be  described  here.  In  addition  the 
Born  and  Rytov  approximations  are  introduced  and  a  linearized  model  for  the 
scattered  field  as  a  function  of  the  object  is  derived. 

In  Chapter  3  the  linearized  wave  equation  is  inverted  to  find  an  expression 
for  the  object  given  the  scattered  field.  This  leads  to  the  Fourier  Diffraction 
Theorem  which  is  fundamental  to  diffraction  tomography.  Finally  several 
experimental  procedures  are  described  that  generate  enough  data  to  uniquely 
determine  the  object. 

Chapter  4  is  a  discussion  of  the  numerical  algorithms  to  invert  the 
scattered  data.  Both  of  these  algorithms  are  computationally  very  expensive 
and  the  algorithm  used  will  depend  on  the  architecture  of  the  available 
computer  resources.  In  addition  some  of  the  signal  processing  issues  will  be 
discussed  and  simulation  results  presented. 
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The  limitations  of  the  first  order  algorithms  are  presented  in  Chapter  5. 
Both  the  mathematical  approximations  and  the  experimental  limitations 
contribute  to  the  error  in  the  final  image  but  in  different  ways.  The 
mathematical  approximations  are  only  valid  for  a  si.. all  range  of  object  and  if 
these  limits  are  exceeded  then  no  amount  of  data  will  improve  the 
reconstruction.  The  experimental  limitations,  on  the  other  hand,  are  entirely 
caused  by  the  ability  to  only  collect  a  finite  amount  of  data.  The  experimental 
errors  can  always  be  reduced  by  using  more  data  or  more  accurate  signal 
processing  algorithms. 

The  severe  limitations  of  first  order  diffraction  algorithms  is  addressed  in 
Chapters  6  and  7.  The  major  problem  in  diffraction  tomography  is  to  find  a 
method  to  invert  the  wave  equation.  In  Chapters  2,  3  and  4  of  this  work  this 
is  done  by  linearizing  the  wave  equation  but  as  seen  this  severly  limits  the 
objects  that  can  be  imaged. 

Chapter  6,  therefore,  discusses  two  approaches  to  model  the  scattered  field 
given  the  (complex)  refractive  index  of  the  object.  This  is  the  forward  problem 
and  both  approaches  are  iterative.  The  simpler  of  the  two  approaches  includes 
more  than  just  the  linear  terms  in  the  perturbation  approach  described  in 
Chapter  2.  This  gives  a  series  solution  for  the  scattered  field  and  simulations 
studying  the  type  of  objects  for  which  these  series  converge  will  be  presented. 
The  second  approach  to  solve  the  forward  problem  exploits  the  simple 
structure  of  the  problem  to  compute  a  brute  force  solution.  Objects  with  large 
refractive  indices  eventually  cause  this  algorithm  to  converge  too  slowly  for  the 
method  to  be  practical. 

Finally  Chapter  7  presents  a  survey  of  several  approaches  that  have  been 
proposed  as  better  solutions  for  the  inverse  problem.  Each  of  these  algorithms 
has  limitations  and  some  of  these  limitations  and  computational  aspects  will  be 
discussed. 
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CHAPTER  2 

DIFFRACTED  PROJECTIONS 


2.1  Introduction 

Tomography  with  diffracting  energy  can  not  he  modeled  with  the  same 
equations  used  to  model  projections  in  conventional,  straight  ray,  tomography. 
Acoustic  and  electromagnetic  waves  do  not  travel  along  straight  rays  and  the 
projections  are  not  line  integrals.  Instead  the  flow  of  energy  will  be  described 
with  the  wave  equation  and  in  the  limit  of  very  short  wavelengths  or  objects 
where  the  effects  of  refraction  are  small  it  will  be  shown  that  the  diffracted 
projections  can  be  approximated  by  a  non  diffracting  projection 

First  consider  the  propagation  of  waves  in  homogeneous  media.  The  wave 
equation  is  a  second  order  linear  differential  equation  and  under  certain 
conditions  it  can  be  shown  that  an  expression  for  the  field  at  every  other  point 
in  space  can  be  written. 

The  problem  is  not  to  image  a  homogeneous  media  but  one  that  is 
inhomogeneous.  To  solve  the  inhomogeneous  wave  equation,  one  of  two 
approximations,  the  Born  or  the  Rytov,  must  be  used.  With  these  two 
approximations  expressions  for  the  field  scattered  by  the  inhomogeneities  of  the 
media  can  be  written. 

The  theory  to  be  discussed  will  be  applicable  to  both  two  and  three 
dimensional  structures.  Even  in  a  three  dimensional  world  a  two  dimensional 
model  can  often  be  used  if  the  object  varies  slowly  in  one  direction.  This 
assumption,  for  example,  is  often  made  in  conventional  computerized 
tomography  where  images  are  made  of  a  single  slice  of  the  object.  The  theory 
of  diffraction  tomography  will  be  presented  almost  entirely  in  two  dimensions 
for  two  reasons.  More  importantly,  the  ideas  behind  the  theory  are  often  easier 
to  visualize  (and  certainly  to  draw)  in  two  dimensions.  In  addition  tec -nology 
has  yet  to  make  it  practical  to  implement  large  three  dimensional  transforms 
and  then  to  display  the  results.  This  limitation  will  certainly  be  eliminated  in 
the  near  future  and  where  the  differences  are  significant  both  the  two  and  three 
dimensional  solutions  will  be  indicated. 
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2.2  Homogeneous  Wave  Equation 

In  a  constant  or  homogeneous  media  the  propagation  of  acoustic  or 
electromagnetic  waves  can  be  modeled  with  the  scalar  Helmholtz  equation.  For 
a  temporal  frequency  of  u>  radians  per  second  (rps)  a  field,  u(7),  satisfies  the 
equation 

(V2  +  k2)u(r)  =0.  (2.1) 

For  homogeneous  media  the  wavenumber,  k0,  is  a  constant  related  to  the 
wavelength.  X,  of  the  wave  by 

2  7T 

k0  =  —  •  (2  2) 

The  wavelength,  X,  is  related  to  the  temporal  frequency  of  the  wave  by  the 
propagation  speed  in  the  media,  c,  or 

2  7T 

X  =  —  c  (2.3) 

u 


Since  the  theory  of  diffraction  tomography  is  normally  derived  based  on 
coherent  fields  the  time  dependence  of  most  fields  will  be  suppressed  in  this 
work.  Thus  all  fields  should  be  multiplied  by  e  Ju"  to  find  the  measured  field  as 
a  function  of  time.  The  extension  of  this  theory  to  broadband  fields  is 
discussed  in  Section  3.4.3 

For  acoustic  (or  ultrasonic)  tomography,  u(7)  can  be  the  pressure  field  at 
position  7.  For  the  electromagnetic  case,  assuming  the  applicability  of  a  scalar 
propagation  equation,  ufr)  may  be  set  equal  to  the  complex  amplitude  of  the 
electric  field  along  its  polarization.  In  both  cases  the  time  dependence  of  the 
fields  are  suppressed  and  u(7)  represents  the  complex  amplitude  of  the  field.  As 
a  function  of  time  and  space  the  field  is  given  by 


u(7,t)  -  Real  Part 


u(7)e“iwl 


(2-4) 


The  vector  gradient  operator,  V,  can  be  expanded  into  its  two 
dimensional  representation  and  the  wave  equation  becomes 

£^i+^a+k2u(T)  =  0,  (2.5) 

ax*  ay* 

As  a  trial  solution  let 

ufr)  =  e^?  (2.6) 

whore  the  vector  k  ~  (kx,ky)  is  the  two  dimensional  propagation  vector  and 
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u(T)  represents  a  two  dimensional  plane  wave  of  spatial  frequency  J  k  |  .  This 
form  of  u(7)  represents  the  basis  function  for  the  two  dimensional  Fourier 
transform;  using  it,  any  two  dimensional  function  can  be  represented  as  as  a 
weighted  sum  of  plane  waves.  Calculating  the  derivatives  as  indicated  in 
equation  (2.5)  it  can  be  seen  that  all  plane  waves  that  satisfy  the  condition 

|7|2  =  k2  +  k2=k02  (2.7) 

are  valid  solutions  to  the  wave  equation.  This  condition  is  consistent  with  an 
intuitive  picture  of  a  wave  and  description  of  the  wave  equation  above,  since 
for  any  frequency  wave  only  a  single  wavelength  can  exist  no  matter  which 
direction  it  propagates. 

The  homogeneous  wave  equation  is  a  linear  differential  equation  so  the 
general  solution  can  be  written  as  a  weighted  sum  of  each  possible  plane  wave 
solution.  In  two  dimensions,  at  a  temporal  frequency  of  u>,  the  field,  u(F)  is 
given  by 

OO  00 

ufr)  =  —  /  o(ky)ej(k*x+k'y)dky  +  /  /?(ky)ej(  M  +  k’y)dky  (2.8) 

00  “OO 

and  by  equation  (2.7) 


The  form  of  this  equation  might  be  surprising  to  the  reader  for  two  reasons. 
First,  the  integral  has  been  split  into  two  parts.  The  coefficients  of  waves 
traveling  to  the  right  are  represented  by  a(ky)  and  those  traveling  to  the  left 
by  /?(ky).  In  addition  the  limits  of  the  integrals  have  been  set  to  go  from  -oo 
to  oo.  For  ky  greater  than  kg  the  radical  in  equation  (2.9)  becomes  imaginary 
and  the  plane  wave  becomes  an  evanescent  wave.  These  are  valid  solutions  to 
the  wave  equation  but  because  ky  is  imaginary  the  exponential  has  a  real  or 
attenuating  component.  This  real  component  causes  the  amplitude  of  the  wave 
to  either  grow  or  decay  exponentially.  In  practice,  these  evanescent  waves  only 
occur  to  satisfy  boundary  conditions,  always  decay  rapidly  far  from  the 
boundary,  and  can  often  be  ignored  at  distance  greater  than  10X  from  the 
inhomogeneity. 

The  limited  range  of  valid  solutions  to  the  wave  equation  allows  (under 
certain  condition)  an  expression  to  be  written  for  the  field  in  all  of  two-space 
given  the  amplitude  of  the  field  along  a  line.  The  three  dimensional  version  of 
this  idea  gives  the  field  in  three-space  if  the  field  is  known  at  all  points  on  a 
plane. 
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Consider  a  source  of  plane  waves  to  the  left  of  a  vertical  line  as  shown  in 
Figure  2.1.  By  calculating  the  one-dimensional  Fourier  transform  of  the  field 
along  the  line  the  field  can  be  decomposed  into  a  number  of  one-dimensional 
components.  Each  of  these  one  dimensional  components  can  then  be  attributed 
to  one  of  the  valid  plane  wave  solutions  to  the  homogeneous  wave  equation 
because  for  any  one  frequency  component,  ky,  there  can  exist  only  two  plane 
waves  that  satisfy  the  wave  equation.  Since  the  incident  field  has  already  been 
constrained  to  propagate  toward  the  right  (all  sources  are  to  the  left  of  the 
measurement  line)  then  a  one-dimensional  Fourier  component  at  a  frequency  of 
ky  can  be  attributed  to  a  two  dimensional  wave  with  a  propagation  vector  of 

(\Ao-k,2,ky) 

This  can  be  put  on  a  more  mathematical  basis  if  the  one-dimensional 
Fourier  transform  of  the  field  is  compared  to  the  general  form  of  the  wave 
equation.  If  waves  that  are  traveling  to  the  left  are  ignored  then  the  general 
solution  to  the  wave  equation  becomes 


OO 

u(C  =  ^/a(kyk-,+l'J,)<lk) 


(2.10) 


Now  if  the  coordinate  system  is  moved  so  that  the  measurement  line  is  at 
x  =  0  then  the  expression  for  the  field  becomes  equal  to  the  one-dimensional 
Fourier  transform  of  the  field  or 


OO 

uffl  =  —  f  a(ky)ejk’ydky. 


(2.11) 


This  equation  establishes  the  link  between  the  one-dimensional  Fourier 
transform  of  the  field  along  the  line  to  the  two-dimensional  field.  The 
coefficients  a(ky)  are  given  from  the  one  dimensional  Fourier  transform  of  the 
field  by 


(2.12) 


a(k  )  -  Fourier  Transform] u(0,y)[. 


The  simple  form  of  a  plane  wave  allows  an  expression  to  be  written 
relating  the  field  on  two  parallel  lines.  If  a  priori  it  is  known  that  all  the 
sources  for  the  field  are  positioned,  for  example,  left  of  the  line  at  x=lo  then 

the  field  u(x=l0,y)  can  be  decomposed  into  its  plane  wave  components.  Given 

if k  I  +  k  v) 

a  plane  wave  *ip|an<?  w»v<-(x  =  Io»y )  =  ,0  the  field  undergoes  a  phase  shift 

as  it  propagates  to  the  line  x—  1(,  and  the  field  can  be  written 


Uplane  wave(x=li,y)  =  aej(k^+kjy)ejMlrW  -  uplane  wave(x  =l0,y)ejkx(trl<))  (2. 13) 

Thus  the  complex  amplitude  of  the  plane  wave  at  x=l,  is  related  to  its 
complex  amplitude  at  x=l0  by  a  factor  of  e^1'  *°). 

The  complete  process  of  finding  the  field  at  a  line  x=l,  follows  in  three 
steps. 

1)  Take  the  Fourier  transform  of  u(x=l<),y)  to  find  the  Fourier 
decomposition  of  u  as  a  function  of  ky. 

2)  Propagate  each  plane  wave  to  the  line  x  =lj  by  multiplying  its 
complex  amplitude  by  the  phase  factor  e^'1  lo^,  where  as  before 


K  =  \Ao~kr 


3)  Find  the  inverse  Fourier  transform  of  the  plane  wave  decomposition 
to  find  the  field  at  u(x=l1,y). 


2.3  Inhomogeneous  Wave  Equation 

For  imaging  in  an  inhomogeneous  media  a  more  general  form  of  the  wave 
equation  is  written  as 


|  V2  +  k(f)2]u(7)  =  0.  (2.14) 

For  the  electromagnetic  case  it  is  necessary  to  ignore  the  effects  of  polarization 
so  that  k(r)  is  a  scalar  function  representing  the  refractive  index  of  the 
medium.  Now  write 

k(f)  =  k0n(f)  =  k0[l  +n^7)]  (2.15) 

where  k0  represents  the  average  wavenumber  of  the  media  and  n^  represents 
the  refractive  index  deviations.  In  general  it  will  be  assumed  that  the  object 
has  a  finite  size  and  therefore  n^(lF)  is  zero  outside  the  object.  Rewriting  the 
wave  equation 

(V2  +  k02)u(T)  =  -k02[n(r)2-l|u(T)  (2.16) 

where  n(T)  is  the  electromagnetic  refractive  index  of  the  media  and  is  given  by 


n(f)  = 


/Wo 


(2.17) 


Here  and  <  have  been  used  to  represent  the  magnetic  permeability  and 
dielectric  constant  and  the  subscript  zero  to  indicate  their  average  values.  This 
new  term,  on  the  right  hand  side  of  equation  (2.16),  is  known  as  a  forcing 
function  for  the  differential  equation  (V2  +  ko)u(r). 


Note  that  equation  (2.16)  is  a  scalar  wave  propagation  equation.  Its  use 
implies  that  there  is  no  depolarization  as  the  electromagnetic  wave  propagates 
through  the  medium.  It  is  known  [Ish78j  that  the  depolarization  effects  can  be 
ignored  only  if  the  wavelength  is  much  smaller  than  the  correlation  size  of  the 
inhomogeneities  in  the  object.  If  this  condition  is  not  satisfied,  then  strictly 
speaking  the  following  vector  wave  propagation  equation  must  be  used 


V2E(F)  +  k02n2E(r)-2V  =  0 


(2.18) 


where  E  is  the  electric  field  vector.  A  vector  theory  for  diffraction  tomography 
based  on  this  equation  has  yet  to  be  developed. 

For  the  acoustic  case,  first  order  approximations  give  the  following  wave 
equation  [Kak84] 

(V2  +  k02)u(7)  =  -k2[n2(f)-l]u(f)  (2.19) 

where  n  is  the  complex  refractive  index  at  position  T,  and  is  equal  to 


n(r)  = 


(2.20) 


where  c0  is  the  propagation  velocity  in  the  medium  in  which  the  object  is 
immersed,  and  c(7)  is  the  propagation  velocity  at  location  ?  in  the  object.  For 
the  acoustic  case  where  only  the  compressional  waves  in  a  viscous  compressible 
fluid  are  involved,  the  speed  of  sound  is  given  by 

c(P)  =  -7=  (2-21) 

V  pK 

where  p  and  k  are  the  local  density  and  the  complex  compressibility  at  location 


The  forcing  function  in  equation  (2.19)  is  only  valid  provided  the  first  and 
higher  order  derivatives  of  the  medium  parameters  can  be  ignored.  If  the 
inhomogeneity  can  be  modeled  as  a  viscous  compressible  fluid,  an  exact  form 
for  the  wave  equation  is  given  by 

(V2  +  k02)u(7)  =  k027Ku  -  V-(Tf.Vu)  (2.22) 


where 


Ik  ~ 


(2.23) 


m 


K 


,« 


'T,  =  y~-  (2  24) 

#c0  and  pQ  are  either  the  compressibility  and  density  of  the  medium  in  which 
the  object  is  immersed,  or  the  average  compressibility  and  the  density  of  the 
object,  depending  upon  how  the  process  of  imaging  is  modeled.  On  the  other 
hand,  if  the  object  is  a  solid  and  can  be  modeled  as  a  linear  isotropic 
viscoelastic  medium,  the  forcing  function  possesses  another  more  complicated 
form.  Since  this  form  involves  tensor  notation,  it  will  not  be  presented  here 
and  the  interested  reader  is  referred  to  (Iwa75). 

Due  to  the  similarities  of  the  electromagnetic  and  acoustic  wave  equations 
a  general  form  of  the  wave  equation  can  be  written  as 

(V2+ko)u(7)  =  -o(7)u(F)  (2.25) 

where 

o(T)  =  koVffH]  (2.26) 

To  hide  some  of  the  mathematical  details  the  term  off)  will  be  used  to 
represent  all  inhomogeneities  of  the  object.  Later  the  object  will  be 
reconstructed  in  terms  of  the  object  function,  off),  and  the  reader  is  referred  to 
equation  (2.26)  to  put  the  reconstruction  in  terms  of  the  refractive  index. 

Consider  the  field,  u(T),  to  be  the  sum  of  two  components.  The  incident 
field,  u0ff),  is  the  field  present  without  any  inhomogeneities  or  a  solution  to  the 
equation 

(V2  +  k02)u0ff)  =  0.  (2.27) 

That  leaves  the  scattered  field,  ugff),  as  that  part  of  the  field  due  to  the  object 
inhomogeneities  or 

us(?)  =  u(f)-u0(f).  (2.28) 

The  wave  equation  becomes 

(V2  +  k02)usff)  =  -u(f)o(f).  (2.29) 

The  scalar  Helmholtz  equation  (2.29)  cannot  be  solved  for  us(f)  directly 
but  a  solution  can  be  written  in  terms  of  the  Green’s  function  (Mor53j.  The 
Green’s  function,  which  is  a  solution  of  the  differential  equation 

(V2  +  k2)gfr|7')  =-6fr-r'), 

is  written  in  three-space  as 


(2.30) 


In  two  dimensions  the  solution  of  (2.30)  is  written  in  terms  of  a  zero-order 
Hankel  function  of  the  first  kind,  and  can  be  expressed  as 

g(F|7')  =  -j-U^koR).  (2.33) 

In  both  cases,  the  Green’s  function,  g(7|7'),  is  only  a  function  of  the  difference 
7~7'  so  the  function  will  often  be  represented  as  simply  g(7-7' ).  Because  the 
object  function  in  equation  (2.30)  represents  a  point  inhomogeneity,  the  Green's 
function  can  be  considered  to  represent  the  field  resulting  from  a  single  point 
scatterer. 

It  is  possible  to  represent  the  forcing  function  of  the  wave  equation  as  an 
array  of  impulses  or 

o(7)u(7)  =  /off'  )u(7'  )6(7-7'  )d7' .  (2.34) 

In  this  equation  the  forcing  function  of  the  inhomogeneous  wave  equation  is 
represented  as  as  a  summation  of  impulses  weighted  by  o(r)u(7)  and  shifted  by 
7.  The  Green’s  function  represents  the  solution  of  the  wave  equation  for  a 
single  delta  function;  because  the  left  hand  side  of  the  wave  equation  is  linear, 
a  solution  can  be  written  by  summing  the  scattered  field  due  to  each  individual 
point  scatterer. 

Using  this  idea,  the  total  field  due  to  the  impulse  o(7'  )u(7'  )6(7-7' )  is 
written  as  a  summation  of  scaled  and  shifted  versions  of  the  impulse  response, 
g(7).  This  is  a  simple  convolution  and  the  total  radiation  from  all  sources  on 
the  right  hand  side  of  (2.29)  must  be  given  by  the  following  superposition: 

us(7)  =  / gfr-7'  )o(7'  )u(7'  )d7'  (2  33) 

At  first  glance  it  might  appear  that  this  is  the  solution  needed  for  the  scattered 
field,  but  it  is  not  that  simple.  An  integral  equation  for  the  scattered  field,  us, 
has  been  written  in  terms  of  the  total  field,  u  =  u0  +  us.  This  equation  needs  to 
be  solved  for  the  scattered  field  and  two  approximations  that  allow  this  to  be 
done  will  now  be  discussed. 


2.4  Approximations  to  the  Wave  Equation 

In  the  last  section  an  inhomogeneous  integral  equation  was  derived  to 
represent  the  scattered  field,  u8(f),  as  a  function  of  the  object,  off).  This 
equation  cannot  be  solved  directly,  but  a  solution  can  be  written  using  either  of 
the  two  approximations  described  here.  These  approximations,  the  Born  and 
the  Rytov,  are  valid  under  different  conditions  but  the  form  of  the  resulting 
solutions  are  quite  similar.  These  approximations  are  the  basis  of  the  Fourier 
Diffraction  Theorem. 

Mathematically  speaking  equation  (2.35)  is  a  Fredholm  equation  of  the 
second  kind.  A  number  of  mathematicians  have  presented  works  describing  the 
solution  of  scattering  integrals  [Hoc73,  Col83]  and  they  should  be  consulted  for 
the  theory  behind  the  approximations  to  be  presented  here. 

2.4.1  The  First  Born  Approximation 

The  first  Born  approximation  is  the  simpler  of  the  two  approaches.  Recall 
that  the  total  field,  uff),  is  expressed  as  the  sum  of  the  incident  field,  u0ff),  and 
a  small  perturbation,  us(T),  or 

uff)  =  u0ff)  +  us(f).  (2.36) 

The  integral  of  equation  (2.35)  is  now  written  as 

usff)  =  Jgfr-7'  )o(T'  )u0(7'  )df'  +/ gff-7'  )o(7'  )usff'  )df'  (2.37) 

but  if  the  scattered  field,  usff),  is  small  compared  to  u0(f)  the  effects  of  the 
second  integral  can  be  ignored  to  arrive  at  the  approximation 

u.,01  “  Wntn  =  JgfM"  )off'  )u0ff'  )df'  (2.38) 

The  first  Born  approximation  is  valid  only  when  the  magnitude  of  the 
scattered  field, 

usff)  =  uff)-u0(f),  (2.39) 

is  smaller  than  the  magnitude  of  the  incident  field,  u0.  If  the  object  is  a 
homogeneous  cylinder  it  is  possible  to  express  this  condition  as  a  function  of 
the  size  of  the  object  and  the  refractive  index.  Let  the  incident  wave,  u0ff),  be 
a  plane  wave  propagating  in  the  direction  of  the  vector,  k0.  For  a  large  object, 
the  field  inside  the  object  will  not  be  well  approximated  by  the  incident  field 

=  UpbWttn  *  AeiC°r  (2.40) 

but  instead  will  be  a  function  of  the  change  in  refractive  index,  n^.  Along  a  ray 


through  the  center  of  the  cylinder  and  parallel  to  the  direction  of  propagation 
of  the  incident  plane  wave  the  field  inside  the  object  becomes  a  slow  (or  fast) 
version  of  the  incident  wave  or 

u„ki,rtm  =  (2.41) 

Since  the  wave  is  propagating  through  the  object  the  phase  difference 
between  the  incident  field  and  the  field  inside  the  object  is  approximately  equal 
to  the  integral  through  the  object  of  the  change  in  refractive  index.  For  a 
homogeneous  cylinder  of  radius  ‘a’  wavelengths  the  total  phase  shift  through 
the  object  becomes 

Phase  Change  =  4^1^  (2.42) 

where  X  is  the  wavelength  of  the  incident  wave.  For  the  Bom  approximation 
to  be  valid,  a  necessary  condition  is  that  the  change  in  phase  between  the 
incident  field  and  the  wave  propagating  through  the  object  be  less  than  7r. 
This  condition  can  be  expressed  mathematically  as  [New66] 

an,<7-  (2  43) 

4 

2.4.2  The  First  Rytov  Approximation 

The  Rytov  approximation  is  another  approximation  to  the  scattered  field 
and  is  valid  under  slightly  different  restrictions.  It  is  derived  by  considering 
the  total  field  to  be  represented  as  a  complex  phase  or  [Ish78] 

(2.4i) 

and  rewriting  the  wave  equation  (2.14) 

(V2  +  k2)u=0  (2.14) 

as 

VV  +  k2e*=0  (2.45) 

V(V^]  +  kV=0  (2.46) 

Vfye*  +  (V0)V  +  kV  -  o  (2.47) 

and  finally 

(V<»2  +  V2<*  +  k02  =  -o(7).  (2.48) 

(Although  all  the  fields  (*'•  and  <t>)  are  a  function  of  7.  to  simplify  the  notation 


^ 
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the  argument  of  these  functions  will  be  dropped.)  Expressing  the  total  complex 
phase,  <p ,  can  be  expressed  as  the  sum  of  the  incident  phase  function,  0O,  and 
the  scattered  complex  phase,  <ps,  or 

=W)  +  «?)  (2.40) 

where 

u0f f)  =  (2.50) 

to  find  that 

(V«>0)2+2V^V^  +  (V^)2  +  VVo  +  VV,+k„2+o(r1  =0.  (2.51) 

As  in  the  Born  approximation  it  is  possible  to  set  the  zero  perturbation 
equation  equal  to  zero.  Doing  this,  the  homogeneous  wave  equation  can  be 
written 

ko +(v<^o)2  +  v2^o  “  (2.52) 

Substituting  this  into  equation  (2.51)  the  wave  equation  becomes 

2V0O  ws  +  VVS  =  -(V<y2-o(7).  (2.53) 

This  equation  is  still  inhomogeneous  but  can  be  linearized  by  considering 
the  following  relation: 

v2(Ms)  =  v(Vuo-^s  +  uoV0s)  (2.54) 

or  by  expanding  the  first  derivative  on  the  right  hand  side  of  this  equation 

V2(uo0s)  =  V2uo-0g  +  2Vuo-V0s  +  uoVVs  (2.55) 

Using  a  plane  wave  for  the  incident  field, 

u0  =  Ae^°r,  (2.56) 

the  second  gradient  of  the  incident  field  is 

V2u0  =  -k02u0  (2.57) 

so  that  equation  (2.55)  may  be  rewritten  as 

2u0V<£0V^s  +  u0VVs  =  V2(u0<Ag)  +  k02u0^s.  (2.58) 

This  result  can  be  substituted  into  equation  (2.53)  to  find 

(V2  +  k2K<.,  =  -u0[(V^)2  +  o(T)]  (2.59) 

The  solution  to  this  differential  equation  can  again  be  expressed  as  an  integral 
equation.  This  becomes 


(2.60) 


Ms  =  /gfF-T')u0[(V^s)2+o(r)jdr'. 

Using  the  Rytov  Approximation  it  is  necessary  to  assume  that  the  term  in 
brackets  in  the  above  equation  can  be  approximated  by 

(Ws)2~o(f)  =*  -off).  (2.61) 

When  this  is  done,  the  first  order  Rytov  approximation  to  the  function  u0<£8 
becomes 


Ms  =  /gf^'Juoff'Jofr'Jdr'.  (2.62) 

v* 

Thus  <ps,  the  complex  phase  of  the  scattered  field,  is  given  by 

*,(?)  =  /  StM"  )"oP"  wr<  )dr* .  (2.63) 

Ml  v 

Substituting  the  expression  for  uB  given  in  equation  (2.38)  the  first  Rytov 
approximation  can  be  written 

m  =  (2.64) 

«o(?) 

The  Rytov  approximation  is  valid  under  a  less  restrictive  set  of  conditions 
than  the  Born  approximation  [Che60,  Kel69).  In  deriving  the  Rytov 
approximation  it  was  necessary  to  assume  that 

/gfP-T'  )u0(F'  )o(f'  )dr*  »  /gfT-r'  )u0(T'  )( WJ2dr' .  (2.65) 

v  v 

If  the  object  is  smaller  then  a  wavelength  then  both  the  field  and  the  object 
can  be  assumed  to  be  constant  compared  to  the  object  function  and  the  above 
relation  can  be  written 

g(TM))u0(0)/o(f')dr'  »  g(TM))u0(O)/(V<£s)2dr'.  (2.66) 

When  the  term  (V<£g)2  is  small  outside  the  object  this  relation  can  be  further 
simplified  to  find 

o(r)  »  (V^)2. 

If  off)  is  written  in  terms  of  the  change  in  refractive  index 


(2.67) 


(2.26) 


o(F)=k02(n2(f)-l]  =k2[(l+n*(r))2-l] 
and  the  square  of  the  refractive  index  is  expanded  to  find 

o(?)  =k02[(l+2ajr)+Ds2(r))-l] 
o(7)  =  k02[2n6(7)  +  n62(?)]. 


(2.68) 

(2.69) 


To  a  first  approximation  the  object  function  is  linearly  related  to  the  refractive 
index  or 


o(7)  »  2k02n*(F). 


(2.70) 


The  condition  needed  for  the  Rytov  approximation  (see  equation  (2.67)  can  be 
rewritten  as 


na  » 


W.)! 


(2.71) 


This  can  be  justified  by  observing  that  to  a  first  approximation  the 
scattered  phase,  <f>3,  is  linearly  dependent  on  the  refractive  index  change,  n$, 
and  therefore  the  first  term  in  equation  (2.65)  above  can  be  safely  ignored  for 
small  n6. 

The  term  V<f>3  is  the  change  in  the  complex  scattered  phase  per  unit 
distance  and  by  dividing  by  the  wavenumber 


a  necessary  condition  for  the  validity  of  the  Rytov  approximation  is 


n*  » 


V<*8X 


(2.72) 


(2.73) 


Unlike  the  Born  approximation,  it  is  the  change  in  scattered  phase,  <f>s,  over  one 
wavelength  that  is  important  and  not  the  total  phase.  Thus,  because  of  the  V 
operator,  the  Rytov  approximation  is  valid  when  the  phase  change  over  a 
single  wavelength  is  small. 

Estimating  uB(F)  for  the  Rytov  case  is  slightly  more  difficult.  In  an 
experiment  the  total  field,  u(T),  is  measured.  An  expression  for  uB(7)  is  found 
by  recalling  the  expression  for  the  Rytov  solution  to  the  total  wave 

u(7)  =  u0  +  us(7)  =  e*0+*• 
and  then  rearranging  the  exponentials  to  find 


(2.74) 
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—  e^o+ 


ug  =  e^le^l 


u»  -  uo(e*'"4 


(2.75) 


(2.76) 


(2.77) 


Inverting  this  to  find  an  estimate  for  the  scattered  phase,  <ps,  the  scattered 
phase  is 


^sOI  =  ln 


us 

— +  1 
u0 


(2.78) 


Then  expand  4>s  in  terms  of  equation  (2.64)  to  obtain  the  following  estimate  for 
the  Rytov  estimate  of  uBfr) 


untn  =  uof1')111 


— +  1 


(2.79) 


Since  the  natural  logarithm  is  a  multiple  valued  function,  one  must  be  careful 
at  each  position  to  choose  the  correct  value.  For  continuous  functions  this  is 
not  difficult  because  only  one  value  will  satisfy  the  continuity  requirement.  On 
the  other  hand  for  discrete  (or  sampled)  signals  the  choice  is  not  nearly  as 
simple  and  one  must  resort  to  a  phase  wrapping  algorithm  to  choose  the  proper 
phase.  Phase  unwrapping  has  been  described  in  a  number  of  works  [Tri77, 
OCo78,  Kav84,  McG82,  Kav84j.  Due  to  the  “  +  1”  factor  inside  the 
logarithmic  term,  this  is  only  a  problem  if  us  is  on  the  order  of  or  larger  than 
u0.  Thus  both  the  Dorn  and  the  Rytov  techniques  can  be  used  to  estimate 

uBfO. 

While  the  Rytov  approximation  is  valid  over  a  larger  class  of  objects,  it  is 
possible  to  show  that  the  Born  and  the  Rytov  approximations  produce  the 
same  result  for  objects  that  are  small  and  deviate  only  slightly  from  the 
average  refractive  index  of  the  medium.  Consider  first  the  Rytov  expression  to 
the  total  field.  This  is  given  by 

u(F)=e*°+*\  (2.80) 

Substituting  an  expression  for  the  scattered  phase,  (2.64)  and  the  incident  field, 
(2.56)  into  this  expression 

u(f)  =  e -jk0sM'+«,Jip(  jk(?r)u,(f)  (2  81) 

or 


>5'v' 

V>-> 


i  iUUULmJ 


u(f)  =  u0ff)e"p(jMr)urf?) 


(2.82) 


For  small  uB,  the  first  exponential  can  be  expanded  in  terms  of  its  pow’r  series. 
Throwing  out  all  but  the  first  two  terms  the  total  field  is  approximately  equal 
to 

u(?)  =  u0(T)[l  +e  jltoS  'ugff)]  (2.83) 

or 

u(7)  =  Uo(7)  +  uBCP).  (2.84) 

Thus  when  the  magnitude  of  the  scattered  field  is  small  the  Rytov  solution  is 
approximately  equal  to  the  Born  solution  given  in  equation  (2.38). 

The  similarity  between  the  expressions  for  the  first  order  Born  and  Rytov 
solutions  will  form  the  basis  of  the  reconstruction  algorithms  to  be  derived 
here.  In  the  Born  approximation  the  complex  amplitude  of  the  scattered  field 
is  measured  and  this  is  used  ^s  an  estimate  of  the  function  uB  while  in  the 
Rytov  case  uB  is  estimated  from  the  complex  phase  of  the  scattered  field. 
Since  the  Rytov  approximation  is  considered  more  accurate  than  the  Born 
approximation  it  should  provide  a  better  estimate  of  uB.  In  Chapter  5  of  this 
work,  after  deriving  reconstruction  algorithms  based  on  the  Fourier  Diffraction 
Theorem,  simulations  comparing  the  Born  and  the  Rytov  approximations  will 
be  discussed. 
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CHAPTER  3 

THE  FOURIER  DIFFRACTION  THEOREM 


3.1  Introduction 

Fundamental  to  diffraction  tomography  is  the  Fourier  Diffraction 
Projection  Theorem,  which  relates  the  Fourier  transform  of  the  measured 
forward  scattered  data  with  the  Fourier  transform  of  the  object.  The  theorem 
is  valid  when  the  inhomogeneities  in  the  object  are  only  weakly  scattering  and 
can  be  stated  as  [Kak84]: 

When  an  object,  f(x,y),  is  illuminated  with  a  plane  wave  as  shown 
in  Figure  3.1,  the  Fourier  transform  of  the  forward  scattered 
fields  measured  on  line  TT’  gives  the  values  of  the  2-D  transform, 
F(cJj,u;2),  of  the  object  along  a  circular  arc  in  the  frequency 
domain,  as  shown  in  the  right  half  of  the  figure. 

The  importance  of  the  theorem  is  made  obvious  by  noting  that  if  an  object  is 
illuminated  by  plane  waves  from  many  directions  over  360°,  the  resulting 
circular  arcs  in  the  (w^u^-plane  will  fill  the  frequency  domain.  The  function 
f(x,y)  may  then  be  recovered  by  Fourier  inversion. 

Before  giving  a  short  proof  of  the  theorem,  first  a  few  words  about  the 
dimensionality  of  the  object  compared  to  that  of  the  scattered  fields.  Although 
the  theorem  talks  about  a  two-dimensional  object,  what  is  actually  meant  is  an 
object  that  does  not  vary  in  the  z  direction.  In  other  words,  the  theorem  is 
about  any  cylindrical  object  whose  cross-sectional  distribution  is  given  by  the 
function  f(x,y).  The  forward  scattered  fields  are  measured  on  a  line  of 
detectors  along  TT’  in  Figure  3.1. 

If  a  truly  three-dimensional  object  is  illuminated  by  a  plane  wave,  the 
forward  scattered  fields  would  now  have  to  be  measured  by  a  planar  array  of 
detectors.  The  Fourier  transform  of  the  fields  measured  by  such  an  array 
would  give  the  values  of  the  3-D  transform  of  the  object  over  a  spherical 
surface.  This  was  first  shown  by  Wolf  [Wol69],  A  more  recent  exposition  is  in 
[Nah84  and  Dev84],  where  the  authors  have  also  presented  a  new  synthetic 
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aperture  procedure  for  a  full  three  dimensional  reconstruction  using  only  two 
rotational  positions  of  the  object.  This  chapter,  however,  will  continue  to  work 
with  two  dimensional  objects  in  the  sense  described  here.  A  recent  work 
describing  some  of  the  errors  in  this  approach  is  [LuZ84j. 


3.2  Decomposing  the  Green’s  Function 

Earlier  in  this  work,  the  scattered  field  due  to  a  weakly  scattering  object 
was  expressed  as  the  convolution 

ub(7)  =  / o{7'  )u0{7 '  )g(F-7'  )dr '  (3. 1 ) 

where  uB(7)  represents  the  complex  amplitude  of  the  field  as  in  the  Born 
approximation  or  the  incident  field,  u0(7),  times  the  complex  scattered  phase, 
^s(7),  in  the  Rytov  approximation.  From  this  integral  there  are  two 
approaches  to  the  derivation  of  the  Fourier  Diffraction  Theorem.  Many 
researchers  [Mue79,  Gre78,  Dev82]  have  expanded  the  Green’s  function  into  its 
plane  wave  decomposition  and  then  noticed  the  similarity  of  the  resulting 
expression  and  the  Fourier  transform  of  the  object.  Alternatively,  if  the 
Fourier  transform  of  each  component  of  this  equation  (3.1)  is  taken  then  the 
Fourier  Diffraction  Theorem  can  be  derived  in  a  manner  that  can  be  easily 
visualized  and  points  towards  efficient  computer  implementations.  This  work 
will  present  both  approaches  to  the  derivation  of  the  Fourier  Diffraction 
Theorem:  the  first  because  the  math  is  more  straightforward,  the  second 
because  it  provides  more  insight  into  the  difference  between  transmission  and 
reflection  tomography. 

First  the  Green’s  function  will  be  decomposed  into  its  plane  wave 
components. 


3.2.1  Plane  Wave  Approach 

The  integral  equation  for  the  scattered  field  (3.1)  can  be  considered  as  a 
convolution  of  the  Green’s  Function,  g(T~T' ),  and  the  product  of  the  object 
function,  o(7' ),  and  the  incident  field,  u0ff '  )■  Consider  the  effect  of  a  single 
plane  wave  illuminating  an  object.  The  forward  scattered  field  will  be 
measured  at  the  receiver  line  as  is  shown  in  Figure  3.2. 

A  single  plane  wave  in  two  dimensions  can  be  represented  as 

u0(7)  =  eJK  r'  (3.2) 

where  R  =  (kx,ky)  satisfies  the  following  relationship 
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k02=k2  +  ky2.  (3.3) 

From  earlier  in  this  work,  the  two  dimensional  Green's  function  is  given 


gtF|T')  =  ^HolkolT-T'l)  (3.4) 

and  H0  is  the  zero-order  Hankel  function  of  the  first  kind.  The  function  H0  has 
the  plane  wave  decomposition  [Mor53] 

CO 

H0(k|7-7'  | )  =  -  /  -ieiN*  *')+^|y  >'Udo  (3  5) 

*■-00? 

where  7  =  (x,y),  7'  =  (x'  ,y' )  and 


0  =  Vk02-a2  (3.6) 

Basically,  equation  (3.5)  expresses  a  cylindrical  wave,  H0,  as  a  superposition  of 
plane  waves.  At  all  points,  the  wave  centered  at  7'  is  traveling  outward;  for 
points  such  that  y>y  '  the  plane  waves  propagate  upward  while  for  y<y  '  the 
plane  waves  propagate  downward.  In  addition,  for  |o|  <k0,  the  plane  waves 
are  of  the  ordinary  type,  propagating  along  the  direction  given  by  tan_,(/?/a). 
However,  for  |  a|  >k0,  0  becomes  imaginary,  the  waves  decay  exponentially 
and  they  are  called  evanescent  waves.  Evanescent  waves  are  usually  of  no 
significance  beyond  about  10  wavelengths  from  the  source. 

Substituting  this  expression,  (3.5),  into  the  expression  for  the  scattered 
field,  (3.1),  the  scattered  field  can  now  be  written 

.  00 

„b(7)  =  -j— /o(r')u0f7')/  “oM*  *')+^  y'Udftdf1  (3.7) 

4k  -oo* 

In  order  to  show  the  first  steps  in  the  proof  of  this  theorem,  assume  for 
notational  convenience  that  the  direction  of  the  incident  plane  wave  is  along 
the  positive  y-axis.  Thus  the  incident  field  is  given  by 

u0fr)  =  ejsV  (3  S) 

where  "sq  =  (O,k0).  Since  in  transmission  imaging  the  scattered  fields  are 
measured  by  a  linear  array  located  at  y  =  10,  where  l0  is  greater  than  any  y- 
coordinate  within  the  object  (see  Figure  3.2),  the  term  j  y— y  '  |  in  the  above 
expression  may  simply  be  replaced  by  l0-y  '  and  the  resulting  form  may  be 
rewritten 
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uB(x,y=lo)  =  “■  /  daf—*~r~^e . "e 


°fjtileiM*  *')  +  ^(io-y  ')lejkay  ' 
0 


(3.9) 


Recognizing  part  of  the  inner  integral  as  the  two-dimensional  Fourier 
transform  of  the  object  function  evaluated  at  a  frequency  of  (o,/?~k0)  the 
scattered  field  can  be  written 


Mx.y-W  =  -j~  /  ^ej(,vx+/?lolO(a,/?-k0Jda 


(3.10) 


where  O  has  been  used  to  designate  the  two  dimensional  Fourier  transform  of 
the  object  function. 

Let  Uq(w,10)  denote  the  Fourier  transform  of  the  one  dimensional  scattered 
field,  U|j(x,l0),  with  respect  to  x,  that  is 


L’bMo)  -  /  un(x,l0)e  Ju,xdx 


(311) 


As  mentioned  before,  the  physics  of  wave  propagation  dictate  that  the  highest 
angular  spatial  frequency  in  the  measured  scattered  field  on  the  line  y=lo  is 
unlikely  to  exceed  k0.  Therefore,  in  almost  all  practical  situations,  Us(u;,l0)  =  0 
for  |  xj  >  k0.  This  is  consistent  with  neglecting  the  evanescent  modes  as 
described  earlier. 

If  the  Fourier  transform  of  the  scattered  field  is  found  by  substituting 
equation  (3.10)  into  equation  (3.11)  then  using  the  following  property  of  Fourier 
integrals 


/  ej(w  °)*dx  =  27r£(u.’-r>) 


(312) 


where  6(  )  is  the  Dirac  delta  function  discussed  in  Chapter  2  the  scattered  field 
can  be  written 


I’pKlo)  = 


± 


2N/k0--o! 


eiv^\0(ft  \/k  2_ft2_ko) 


(3.13) 


for  |  o  |  <  k0. 

This  expression  relates  the  two  dimensional  Fourier  transform  of  the  object  to 
the  one  dimensional  Fourier  transform  of  the  field  at  the  receiver  line.  The 
factor 
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Jv/ko  a4lo 


'k$~a2 


is  a  simple  constant  for  a  fixed  receiver  line.  As  «  varies  from  ~k0  to  k0,  the 
coordinates  (a,  >A  0-a2-k0)  in  the  Fourier  transform  of  the  object  function 
trace  out  a  semicircular  arc  in  the  (u,v)-plane  as  shown  in  Figure  3.1.  This 
proves  the  Fourier  Diffraction  Theorem. 

To  summarize,  if  the  Fourier  transform  of  the  forward  scattered  data  is 
found  when  the  incident  illumination  is  propagating  along  the  positive  y-axis, 
the  resulting  transform  will  be  zero  for  angular  spatial  frequencies  |  or |  >k0 
For  J  or |  < k0,  the  transform  of  the  data  gives  values  of  the  Fourier  transform 
of  the  object  on  the  semicircular  arc  are  shown  in  Figure  3.1  in  the  (u,v)-plane. 
The  endpoints  A  and  B  of  the  semicircular  arc  are  at  a  distance  of  v/2k0  from 
the  origin  in  the  frequency  domain. 

3.2.2  Fourier  Transform  Approach 

Another  approach  to  the  derivation  of  the  Fourier  Diffraction  Theorem  is 
possible  if  the  scattered  field 

uBfr)  =  / of7'  )uo(7'  tef7-7'  W  (315) 

is  considered  entirely  in  the  Fourier  domain.  The  plots  of  Figure  3.3  will  be 
used  to  illustrate  the  various  transformations  that  take  place. 

Again  consider  the  effect  of  a  single  plane  wave  illuminating  an  object. 
The  forward  scattered  field  will  be  measured  at  the  receiver  line  as  is  shown  in 
Figure  3.2. 

The  integral  equation  for  the  scattered  field  can  be  considered  as  a 
convolution  of  the  Green’s  Function,  g(T-f' ),  and  the  product  of  the  object 
function,  off'),  and  the  incident  field,  u0(f ' ).  First  define  the  following  Fourier 
transform  pairs. 

off)  ~  Q(i\) 


gfr-f')-G(R) 


uff)  «-  U(l\) 


(3.16) 


The  integral  solution  to  the  wave  equation  can  now  be  written  in  terms  of 
these  Fourier  transforms  or 


(3.17) 


Ua(A)  =  G(A)jO(A)*U0(A)J 

where  **’  has  been  used  to  represent  convolution  and  A  —  (0,7).  In  equation 
(3.2)  an  expression  for  u0  was  presented.  It’s  Fourier  transform  is  given  by 


U0(A)  =  2jt«(A-K) 


(3.18) 


and  thus  the  convolution  of  equation  (3.17)  becomes  a  shift  in  the  frequency 
domain  or 

O(A)*U0(f)  =  2ttO(A-R).  (3.19) 

This  convolution  is  illustrated  in  Figures  3.3a-c  for  a  plane  wave  propagating 
with  direction  vector,  R  =  (0,k0).  Figure  3.3a  shows  the  Fourier  transform  of  a 
single  cylinder  of  radius  IX  and  Figure  3.3b  is  the  Fourier  transform  of  the 
incident  field.  The  resulting  multiplication  in  the  space  domain  or  convolution 
in  the  frequency  domain  is  shown  in  Figure  3.3c. 

To  find  the  Fourier  transform  of  the  Green’s  function  the  Fourier 
transform  of  the  equation  for  a  point  scatterer 

(V2  +  k02)gfr|-r' )  =  -6(r-r  ),  (3.20) 


is  taken  to  find 


(-A2+k02)G(A|7')  =  -e~i*?'. 


(3.21) 


Rearranging  terms  the  following  expression  for  the  Fourier  transform  of  the 
Green’s  function  is  found 


G(A|T')  = 


A2-kn 


This  has  a  singularity  for  all  A  such  that 

|A|2=o2  +  72  =  k02 


(3.22) 


(3.23) 


An  approximation  to  G(A)  is  shown  in  Figure  3.3d. 

The  Fourier  transform  representation  is  misleading  because  it  represents  a 
point  scatterer  as  both  a  sink  and  a  source  of  waves.  A  single  plane  wave 
propagating  from  left  to  right  can  be  considered  in  two  different  ways 
depending  on  the  point  of  view.  From  the  left  side  of  the  scatter,  the  point 
scatterer  represents  a  sink  to  the  wave  while  to  the  right  of  the  scatterer  the 


wave  is  spreading  from  a  source  point.  Clearly,  it  is  not  possible  for  a  scatterer 
to  be  both  a  point  source  and  sink,  and  later  when  our  expression  for  the 
scattered  field  is  inverted,  it  will  be  necessary  to  choose  a  solution  that  leads  to 
outgoing  waves  only. 

The  effect  of  the  convolution  shown  in  equation  (3.15)  is  a  multiplication 
in  the  frequency  domain  of  the  shifted  object  function,  (3.19),  and  the  Green’s 
function,  (3.22),  evaluated  at  7'  =  0.  The  scattered  field  is  written  as 

u„(A>  =  <3.24) 

This  result  is  shown  in  Figure  3.3e  for  a  plane  wave  propagating  along  the  y- 
axis.  Since  the  largest  frequency  domain  components  of  the  Green’s  function 
satisfy  equation  (3.23),  the  Fourier  transform  of  the  scattered  field  is  dominated 
by  a  shifted  and  sampled  version  of  the  object’s  Fourier  transform. 

An  expression  for  the  field  at  the  receiver  line  will  now  be  derived.  For 
simplicity  it  will  continue  to  be  assumed  that  the  incident  field  is  propagating 
along  the  positive  y  axis  or  R  =  (0,k0).  The  scattered  field  along  the  receiver 
line  (x,y  =  lo)  is  simply  the  inverse  Fourier  transform  of  the  field  in  equation 
(3.24).  This  is  written  as 

OO  00 

u(x,y=l0)  =  /  J  l)s(A)eiA?dad7  (3.25) 
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which,  using  (3.24),  can  be  expressed  as 


4 ir  -po-oo  c^  +  r-ko 


(3.26) 


First  find  the  integral  with  respect  to  7.  For  a  given  a,  the  integral  has  a 
singularity  for 

7I2  =  ±s/k$-a2  (3.27) 

Using  contour  integration  the  integral  can  be  evaluated  with  respect  to  7  along 

the  path  shown  in  Figure  3.4.  By  adding  — —  of  the  residue  at  each  pole  the 

2  7T 

scattered  field  is  expressed 
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r,  = 


_  jO(a,  y  k^-g^kp)  ]y/^]o 


2\/k  ft2-a2 


(3.29) 


r,  = 


_  -jO(a,-\/ko2-a2-ko)  -jv^^ 


\Fa? 


(3.30) 


Examining  the  above  pair  of  equations  it  can  be  seen  that  represents  the 
solution  in  terms  of  plane  waves  traveling  along  the  positive  y  axis  while  T2 
represents  plane  waves  traveling  in  the  -y  direction. 

As  was  discussed  earlier,  the  Fourier  transform  of  the  Green’s  function 
(3.22)  represents  the  field  due  to  both  a  point  source  and  a  point  sink  but  the 
two  solutions  are  distinct  for  receiver  lines  that  are  outside  the  extant  of  the 
object.  First  consider  the  scattered  field  along  the  line  y  =  lo  where  ^  is 
greater  than  the  y-coordinate  of  all  points  in  the  object.  Since  all  scattered 
fields  originate  in  the  object,  plane  waves  propagating  along  the  positive  y  axis 
represent  outgoing  waves  while  waves  propagating  along  the  negative  y  axis 
represent  waves  due  to  a  point  sink.  Thus  for  y>object  (i.e.  the  receiver  line  is 
above  the  object)  the  outgoing  scattered  waves  are  represented  by  T,  or 


u.lx.yl  =  —;/r,(n;y|ei“"do  y>object 


(3.31) 


Conversely  for  a  receiver  along  a  line  y  =  ^  where  ^  is  less  than  the  y- 
coordinate  of  any  point  in  the  object  the  scattered  field  is  represented  by  f2  or 


us(x,y)  =  “/F2(a;y)ej°xdft  y<object 

Z7T 


(3.32) 


In  general  the  scattered  field  will  be  written  as  as 

«s(x,y)  =  ~-/r(a;y)e*axda  (3.33) 

Z7T 

and  it  will  be  understood  that  values  of  the  square  root  in  the  expression  for  T 
should  be  chosen  that  lead  only  to  outgoing  waves. 

Taking  the  Fourier  transform  of  both  sides  of  equation  (3.33)  the  Fourier 
transform  of  the  scattered  field  at  the  receiver  line  is  written 

/u(x,y=l0)e  JO,xdx  =  F(a,l0).  (3.34) 

But  since  by  equations  (3.29)  and  (3.30),  F(q,10)  is  equal  to  a  phase  shifted 
version  of  the  object  function  then  the  Fourier  transform  of  the  scattered  field 
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along  the  line  y  =  l0  is  related  to  the  Fourier  transform  of  the  object  along  a 
circular  arc.  The  use  of  the  contour  integration  is  further  justified  by  noting 
that  only  those  waves  that  satisfy  the  relationship 

a2  +  72=k02  (3.35) 

will  be  propagated  and  thus  it  is  safe  to  ignore  all  waves  not  on  the  k0-circle. 

This  result  is  diagramed  in  Figure  3.5.  The  circular  arc  represents  the 
locus  of  all  points  (0,7)  such  that  7  =  ±\Jk$-a2.  The  solid  line  shows  the 
outgoing  waves  for  a  receiver  line  at  y  =  l0  above  the  object.  This  can  be 
considered  transmission  tomography.  Conversely  the  dashed  line  indicates  the 
locus  of  solutions  for  the  reflection  tomography  case,  or  y=l0  is  below  the 
object. 

3.3  Limit  of  the  Fourier  Diffraction  Theorem 

While  at  first  the  derivations  of  the  Fourier  Slice  Theorem  and  the  Fourier 
Diffraction  Theorem  seem  quite  different,  it  is  interesting  to  note  that  in  the 
limit  of  very  high  energy  waves  or,  equivalently,  very  short  wavelengths  the 
Fourier  Diffraction  Theorem  is  closely  approximated  by  the  Fourier  Slice 
Theorem.  Recall  that  the  Fourier  transform  of  a  diffracted  projection 
corresponds  to  samples  of  the  two  dimensional  Fourier  transform  of  an  object 
along  a  circular  arc.  As  shown  in  Figure  3.1  the  radius  of  the  arc  is  equal  to  k0 
which  is  given  by 

k0  =  Y  (3.36) 

and  X  is  the  wavelength  of  the  energy.  As  the  wavelength  is  decreased,  the 
wavenumber,  k0,  and  the  radius  of  the  arc  in  the  object’s  Fourier  domain 
grows.  This  process  is  illustrated  in  Figure  3.6  where  the  semicircular  arc 
resulting  from  a  diffraction  experiment  is  shown  at  six  different  frequencies. 

An  example  might  make  this  idea  clearer.  Compare  an  ultrasonic 
diffraction  apparatus  and  a  typical  x-ray  scanner.  The  ultrasonic  experiment 
might  be  carried  out  at  a  frequency  of  5  MHz  and  a  wavelength  in  water  of  .3 
mm.  This  corresponds  to  a  k0  of  333  radians/meter.  On  the  other  hand,  an 
x-ray  source  with  a  100  keV  beam  has  a  wavelength  of  .012  /xM.  The  result  is 
that  a  diffraction  experiment  gives  samples  along  an  arc  of  radius  5xl08 
radians/meter.  Certainly  for  all  physiological  features  (i.e.  resolutions  of  < 
1000  radiar.s/meter)  the  arc  can  be  considered  a  straight  line  and  the  Fourier 
Slice  Theorem  is  an  excellent  model  of  the  propagation  of  x-rays. 
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Estimate  of  the  two  dimensional  Fourier  transform  of  the 
object  are  available  along  the  solid  arc  for  transmission 
tomography  and  the  dashed  arc  for  reflection  tomography. 
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Figure  3.6  As  the  illuminating  frequency  is  increased  the  Fourier 

Diffraction  Theorem  becomes  equivalent  to  the  Fourier  Slice 
Theorem. 
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3.4  The  Data  Collection  Process 


The  best  that  can  be  hoped  for  in  any  tomographic  experiment  is  to 
estimate  the  Fourier  transform  of  the  object  for  all  frequencies  within  a  disk 
centered  at  the  origin.  For  objects  that  do  not  have  any  frequency  content 
outside  the  disc  then  the  reconstruction  procedure  is  perfect. 

There  are  several  different  procedures  that  can  be  used  to  estimate  the 
object  function  from  the  scattered  field.  A  single  plane  wave  provides  exact 
information  (up  to  a  frequency  of  v/2k0)  about  the  Fourier  transform  of  the 
object  along  a  circular  arc.  Two  of  the  simplest  procedures  involve  changing 
the  orientation  and  frequency  of  the  incident  plane  waves  to  move  the 
frequency  domain  arcs  to  a  new  position.  By  appropriately  choosing  an 
orientation  and  a  frequency  it  is  possible  to  estimate  the  Fourier  transform  of 
the  object  at  any  given  frequency.  In,  addition  it  is  possible  to  change  the 
radius  of  the  semicircular  arc  by  varying  the  frequency  of  the  incident  field  and 
thus  generating  an  estimate  of  the  entire  Fourier  transform  of  the  object. 

3.4.1  Plane  Wave  Illumination 

The  most  straightforward  data  collection  procedure  consists  of  rotating  the 
object  and  measuring  the  scattered  field  for  different  orientations.  Each 
orientation  will  produce  an  estimate  of  the  object’s  Fourier  transform  along  a 
circular  arc  and  these  arcs  will  rotate  as  the  object  is  rotated.  When  the  object 
is  rotated  through  a  full  360  degrees  an  estimate  of  the  object  will  be  available 
for  the  entire  Fourier  disk. 

The  coverage  for  this  method  is  shown  in  Figure  3.7  for  a  simple 
experiment  with  8  projections  of  9  samples  each.  Notice  that  there  are  two 
arcs  that  pass  through  each  point  of  Fourier  space.  Generally  it  will  be 
necessary  to  choose  one  estimate  as  better. 

On  the  other  hand  if  the  reflected  data  is  collected  by  measuring  the  field 
on  the  same  side  of  the  object  as  the  source  then  estimates  of  the  object  are 
available  for  frequencies  greater  than  v/2k0.  This  follows  from  Figure  3.5. 

The  first  experimental  results  for  diffraction  tomography  were  presented 
by  Carter  and  Ho  [Car70,  Car74,  Car76  and  Ho76].  They  used  an  optical  plane 
wave  to  illuminate  a  small  glass  object  and  were  able  to  measure  the  scattered 
fields  using  a  hologram.  Later  a  group  of  researchers  at  the  University  of 
Minnesota  carried  out  the  same  experiments  using  ultrasound  and  gelatine 
phantoms.  Their  results  are  discussed  in  [Kav82]. 
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3.4.2  Synthetic  Aperture 

Nahamoo  and  Kak  [Nah82,  Nah84]  and  Devaney  [Dev84j  have  proposed  a 
method  that  requires  only  two  rotational  views  of  an  object.  Consider  an 
arbitrary  source  of  waves  in  the  transmitter  plane  as  shown  in  Figure  3.8.  The 
transmitted  field,  ut,  can  be  represented  as  a  weighted  set  of  plane  waves  by 
taking  the  Fourier  transform  of  the  transmitter  aperture  function  [G0068]. 
Doing  this  the  transmitted  field  can  be  expressed  as 

OO 

«t(x)  =  TT  /  At(kx)ejk‘xdkx.  (3.37) 

4ir  — oo 

Moving  the  source  to  a  new  position,  rj,  the  plane  wave  decomposition  of  the 
transmitted  field  becomes 

ut(x;,f)  =  ~J  f  (a^x )ejk,f,)eJ,CxXdkx.  (3.38) 

47l  -oq 


Given  the  plane  wave  decomposition,  the  incident  field  in  the  plane  follows 
simply  as 


Ui(»?;*>y)  =  / 


4^ 


At(kx)e' 


ej(k,x+kry)dk^ 


(3.39) 


In  equation  (3.34)  an  equation  for  the  scattered  field  from  a  single  plane 
wave  was  presented.  Because  of  the  linearity  of  the  Fourier  transform,  the 
effect  of  each  plane  wave,  e^k,x+lc,y^  can  be  weighted  by  the  expression  in 
brackets  above  and  superimposed  to  find  the  Fourier  transform  of  the  total 
scattered  field  due  to  the  incident  field  ut(x;»j)  as  [Nah82] 

°r  f  it  n)  O(o— kx,7-ky) 

U,(„;o)  =  J  [Al(k,yk-',J - ^dk,.  (3  40) 

Taking  the  Fourier  transform  of  both  sides  with  respect  to  the  transmitter 
position,  t),  the  Fourier  transform  of  the  scattered  field  with  respect  to  both  the 
transmitter  and  the  receiver  position  is  given  by 

0(or— kx,7— kv) 

U,(kx;o)  =  At(k„)  y'k,.  (3  41) 

This  approach  gets  the  name  synthetic  aperture  because  a  phase  is  added 
to  the  field  measured  for  each  transmitter  position  to  synthesize  a  transmitted 
plane  wave.  Thus  this  method  has  a  lot  in  common  with  the  theory  of  phased 
arrays.  Figure  3.9  shows  that  by  properly  phasing  the  wave  transmitted  at 
each  transmitter  location  a  plane  wave  can  be  generated  that  travels  in  an 
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Figure  3.9  By  adding  a  phase  to  the  field  transmitted  from  each 

transmitter  any  desired  plane  wave  can  be  synthesized. 


arbitrary  direction.  Since  the  system  is  linear  it  doesn’t  matter  whether  the 
phase  is  added  to  the  transmitted  signal  or  as  part  of  the  reconstruction 
procedure.  Thus  multiplying  the  received  field  for  each  transmitter  position  by 
the  pure  phase  term  e^kx,,,  where  r)  represents  the  location  of  the  transmitter,  is 
equivalent  to  an  experiment  with  an  incident  plane  wave  with  the  direction 
vector  (kx,^/kQ-kx2). 

By  collecting  the  scattered  field  along  the  receiver  line  as  a  function  of 
transmitter  position,  17,  an  expression  can  be  written  for  the  scattered  field. 
Like  the  simpler  case  with  plane  wave  incidence,  the  scattered  field  is  related  to 
the  Fourier  transform  of  the  object  along  an  arc.  Unlike  the  previous  case, 
though,  the  coverage  due  to  a  single  view  of  the  object  is  a  pair  of  circular 
disks  as  shown  in  Figure  3.10.  Here  a  single  view  consists  of  transmitting  from 
all  positions  in  a  line  and  measuring  the  scattered  field  at  all  positions  along 
the  receiver  line.  By  rotating  the  object  by  90  degrees  it  is  possible  to  generate 
the  complimentary  disk  and  to  fill  the  Fourier  domain. 

The  coverage  shown  in  Figure  3.10  is  constructed  by  calculating  (Iv — A)  for 
all  vectors  (R)  and  (A)  that  satisfy  the  experimental  constraints.  Not  only 
must  each  vector  satisfy  the  wave  equation  but  it  is  also  necessary  that  only 
forward  traveling  plane  waves  be  used.  The  dashed  line  in  Figure  3.10  shows 
the  valid  propagation  vectors  (-A)  for  the  transmitted  waves.  To  each  possible 
vector  (~A)  a  semicircular  set  of  vectors  representing  each  possible  received 
wave  can  be  added.  The  locus  of  received  plane  waves  is  shown  as  a  solid 
semi-circle  centered  at  each  of  the  transmitted  waves  indicated  by  an  ‘x’.  The 
entire  coverage  for  the  synthetic  aperture  approach  is  shown  as  the  shaded 
areas. 

In  addition  to  the  diffraction  tomography  configurations  proposed  by 
Mueller  and  Nahamoo  other  approaches  have  been  proposed.  In  Vertical 
Seismic  Profiling  (VSP)  [Dev84]  the  scattering  between  the  surface  of  the  Earth 
and  a  borehole  is  measured.  Alternately  a  broadband  incident  field  can  be  used 
to  illuminate  the  object.  In  both  cases,  the  goal  is  to  estimate  the  Fourier 
transform  of  the  object. 

In  geophysical  imaging  it  is  not  possible  to  generate  or  receive  waves  from 
all  positions  around  the  object.  If  it  is  possible  to  drill  a  borehole  then  it  is 
possible  to  perform  VSP  and  obtain  information  about  most  of  the  object.  A 
typical  experiment  is  shown  in  Figure  3.11.  So  as  to  not  damage  the  borehole, 
acoustic  waves  are  generated  at  the  surface  using  acoustic  detonators  or  other 
methods  and  the  scattered  field  is  measured  in  the  borehole. 


The  coverage  in  the  frequency  domain  is  similar  to  the  synthetic  aperture 
approach.  Plane  waves  at  an  arbitrary  downward  direction  are  synthesized  by 
appropriately  phasing  the  transmitting  transducers.  The  receivers  will  receive 
any  waves  traveling  to  the  right.  The  resulting  coverage  for  this  method  is 
shown  in  Figure  3.12a.  If  can  be  assumed  that  the  object  function  is  real 
valued  then  the  symmetry  of  Fourier  transform  for  real  valued  functions  can  be 
used  to  obtain  the  coverage  in  Figure  3.12b. 

3.4.3  Broadband  Illumination 

It  is  also  possible  to  perform  an  experiment  for  broadband  illumination 
[Ken82],  Up  until  this  point  only  narrow  band  illumination  has  been 
considered;  wherein  the  field  at  each  point  can  be  completely  described  by  its 
complex  amplitude. 

Now  consider  a  transducer  that  illuminates  an  object  with  a  wave  of  the 
form  at(kk,t).  Taking  the  Fourier  transform  in  the  time  domain  this  wave  can 
be  decomposed  into  a  number  of  experiments.  Let 

OO 

At(kx,w)  =  /  at(kx,t)e“ju'tdt  (3.42) 

-oo 

where  u  is  related  to  kw  by 

(3.43) 

c  is  the  speed  of  propagation  in  the  media  and  the  wavevector  (kx,ky)  satisfies 
the  wave  equation 

kx2  +  ky2=kj-  (3.44) 

If  a  plane  wave  illumination  of  spatial  frequency  kx  and  a  temporal 
frequency  w  leads  to  the  scattered  field  us(kx,w;y)  then  the  total  scattered  field 
is  given  by  a  weighted  superposition  of  the  scattered  fields  or 

OO 

u3(k*;y)  =  /  At(kx,w)us(kx,a;;y)dw.  (3.45) 

-00 

For  plane  wave  incidence  the  coverage  for  this  method  is  shown  in  Figure 
3.13a.  Figure  3.13b  shows  that  by  doing  four  experiments  at  0,  90,  180  and 
270  degrees  it  is  possible  to  gather  information  about  the  entire  object. 
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CHAPTER  4 

RECONSTRUCTION  PROCEDURES 


4.1  Introduction 

The  Fourier  Diffraction  Theorem  as  derived  in  Chapter  3  shows  that  when 
an  object  is  illuminated  with  a  plane  wave  traveling  in  the  positive  y-direction, 
the  Fourier  transform  of  the  forward  scattered  fields  gives  values  of  the  object’s 
Fourier  transform  on  an  arc.  Therefore,  if  an  object  is  illuminated  from  many 
different  directions  it  is  possible,  in  principle,  to  fill  up  a  disc  of  diameter  \/2k 
in  the  frequency  domain  with  samples  of  the  Fourier  transform  of  the  object 
and  then  reconstruct  the  object  by  direct  Fourier  inversion.  Therefore, 
diffraction  tomography,  using  forward  scattered  data  only,  determines  the 
object  up  to  a  maximum  angular  spatial  frequency  of  \/2k.  To  this  extent,  the 
reconstructed  object  is  a  low  pass  version  of  the  original.  In  practice,  the  loss  of 
resolution  caused  by  this  bandlimiting  is  negligible,  being  more  influenced  by 
considerations  such  as  the  aperture  sizes  of  the  transmitting  and  receiving 
elements,  etc. 

The  fact  that  the  frequency  domain  samples  are  available  over  circular 
arcs,  whereas  for  convenient  display  it  is  desired  to  have  samples  over  a 
rectangular  lattice,  is  a  source  of  computational  difficulty  in  reconstruct^  • 
algorithms  for  diffracting  tomography.  It  should  also  be  clear  that  by 
illuminating  the  object  over  360  ° ,  a  double  coverage  of  the  frequency  domain  is 
generated;  note,  however,  that  this  double  coverage  is  uniform.  If  the 
illumination  is  restricted  to  a  portion  of  360°  there  still  will  be  a  complete 
coverage  of  the  frequency  domain,  however  in  that  case  there  would  be  patches 
in  the  (u>i,w2)-plane  where  there  would  be  a  double  coverage.  In  reconstructing 
from  circular  arc  grids  to  rectangular  grids,  it  is  often  easier  to  contend  with  a 
uniform  double  coverage,  as  opposed  to  a  coverage  that  is  single  in  most  areas 
and  double  in  patches. 

However,  for  some  applications  not  given  to  data  collection  from  all 
possible  directions,  it  is  useful  to  bear  in  mind  that  it  is  not  necessary  to  go 
completely  around  an  object  to  get  complete  coverage  of  the  frequency  domain. 
In  principle,  it  should  be  possible  to  get  an  equal  quality  reconstruction  when 
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illumination  angles  are  restricted  to  a  180  plus  interval,  the  angles  in  excess 
of  180  °  being  required  to  complete  the  coverage  of  the  frequency  domain. 

There  are  two  computational  strategies  for  reconstructing  the  object  given 
measurements  of  the  scattered  field.  As  pointed  out  by  [Sou84aj  the  two 
algorithms  can  be  considered  as  interpolation  in  the  frequency  domain  and  in 
the  space  domain  and  are  analogous  to  the  direct  Fourier  inversion  and 
backprojection  algorithms  of  conventional  tomography.  Unlike  conventional 
tomography,  where  backprojection  is  the  preferred  approach,  the 
computational  expense  of  space  domain  interpolation  of  diffracted  projections 
makes  frequency  domain  interpolation  the  preferred  approach. 

The  remainder  of  this  section  will  consist  of  derivations  of  the  frequency 
domain  and  space  domain  interpolation  algorithms.  In  both  cases  plane  wave 
illumination  will  be  assumed  and  the  reader  is  referred  to  [Dev82,  Pan83]  for 
reconstruction  algorithms  for  the  synthetic  aperture  approach  and  to  [Sou84b] 
for  the  general  case. 

4.2  Frequency  Domain  Interpolation 

In  order  to  discuss  the  frequency  domain  interpolation  between  a  circular 
grid  on  which  the  data  is  generated  by  diffraction  tomography,  and  a 
rectangular  grid  suitable  for  image  reconstruction,  parameters  for  representing 
each  grid  must  be  selected  and  then  the  relationship  between  the  two  sets  of 
parameters  written. 

In  Chapter  3,  UB(w)  was  used  to  denote  the  Fourier  transform  of  the 
transmitted  data  when  an  object  is  illuminated  with  a  plane  wave  traveling 
along  the  positive  y  direction.  Now  UB^(u;)  is  used  to  denote  this  Fourier 
transform,  where  the  subscript  <p  indicates  the  angle  of  illumination.  This  angle 
is  measured  as  shown  in  Figure  4.1.  Similarly,  Q {u,<p)  will  be  used  to  indicate 
the  values  of  ofw,,^)  along  a  semi-circular  arc  oriented  at  an  angle  <f>  as  shown 
in  Figure  4.2  or 


2 - 2_l 


w  <k0. 


Therefore,  when  an  illuminating  plane  wave  is  incident  at  angle  <f>,  the  equality 
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can  be  rewritten  as 
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Figure  4.2  A  second  change  of  variables  is  used  to  relate  the  projection 

data  to  the  object’s  Fourier  transform. 
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In  most  cases  the  transmitted  data  will  be  uniformly  sampled  in  space, 
and  a  discrete  Fourier  transform  of  this  data  will  generate  uniformly  spaced 
samples  of  UB^(w)  in  the  u>  domain.  Since  Q(u;)  is  the  Fourier  transform  of  the 
object  along  the  circular  arc  AOB  in  Figure  4.1  and  since  k  is  the  projection  of 
a  point  on  the  circular  arc  on  the  tangent  line  CD,  the  uniform  samples  of  Q  in 
k  translate  into  non-uniform  samples  along  the  arc  AOB  as  shown  in  Figure 
4.3.  For  this  reason  designate  each  point  on  the  arc  AOB  by  its  ( u,<f> ) 
parameters.  [Note  that  (u,<f>)  are  not  the  polar  coordinates  of  a  point  on  arc 
AOB  in  Figure  4.2.  Therefore,  u  is  not  the  radial  distance  in  the  (w,,o;2)  plane. 
For  point  E  shown,  the  parameter  u>  is  obtained  by  projecting  E  onto  line  CD.] 
The  rectangular  coordinates  in  the  frequency  domain  will  remain  (u>,,u;2). 

Before  the  relationships  between  and  (w,,a>2),  is  presented  it  must  be 

mentioned  that  the  points  generated  by  the  AO  and  OB  portions  of  the  arc 
AOB  must  be  considered  separately  as  0  is  varied  from  0  to  27r.  This  is  done 
because  as  mentioned  before,  the  arc  AOB  generates  a  double  coverage  of  the 
frequency  domain,  as  <f>  is  varied  from  0  to  2tt,  which  is  undesirable  for 
discussing  a  one-to-one  transformation  between  the  (u,0)  parameters  and  the 
(u>,,w2)  coordinates. 

Now  reserve  (oj,<I>)  parameters  to  denote  the  arc  grid  generated  by  one 
projection.  It  is  important  to  note  that  for  this  arc  grid,  uj  varies  from  0  to  k 
and  </>  from  0  to  2n. 

The  transformation  equations  between  and  (w  1,oj2)  will  now  be 

presented.  This  is  accomplished  in  a  slightly  round-about  manner  by  first 
defining  polar  coordinates  (0,0)  in  the  (Wj,u>2)-plane  as  shown  in  Figure  4.2.  In 
order  to  go  from  (wj,w2)  to  (u,<f>),  first  transform  from  the  former  coordinates  to 
(0,0)  and  then  from  (0,0)  to  (u>,d>).  The  rectangular  coordinates  (w,,w2)  are 
related  to  the  polar  coordinates  (0,0)  by  (Figure  4.2) 

0  —  \Ju\  +  <^2  (4.4) 


0  -  tan 
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(4.5) 


In  order  to  relate  (0,0)  to  (u\0),  a  new  angle  /?,  which  is  the  angular  position  of 
a  point  (uq,w2)  on  arc  OB  in  Figure  4.2,  is  introduced.  Note  from  the  figure 
that  the  point  characterized  by  angle  0  is  also  characterized  by  parameter  u>. 


Figure  4.3 


Uniformly  sampling  the  projection  in  the  space  domain  leads 
to  uneven  spacing  of  the  samples  of  the  Fourier  transform  of 
the  object  along  the  semi-circular  arc. 


The  relationship  between  w  and  p  is  given  by 

w  =  k  sin/?.  (4.4) 

The  following  relationship  exists  between  the  polar  coordinates  (Q,0)  on  the  one 
hand  and  the  parameters  /?  and  $  on  the  other: 

f  -  (4  6) 


*  =  «  +  f  +  {•  (4.7) 

By  substituting  Equation  (4.6)  in  (4.4)  and  then  using  (4.4),  u  can  be  expressed 
in  terms  of  ul  and  w2.  This  result  is  shown  below.  Similarly,  by  substituting 
Equation  (4.5)  in  (4.7),  the  following  expression  is  obtained  for  u  and  <f> 


_  .  .  JvV  +  W22 

w  -  sin  2  sin  1  - : - 

2k  JJ 


* =  tan" +  “'l  arn +  f  ■ 

These  are  the  transformation  equations  for  interpolating  from  the  (w,0) 
parameters  used  for  data  representation  to  the  (w,,w2)  parameters  needed  for 
inverse  transformation. 

To  convert  a  particular  rectangular  point  into  (w,0)  domain,  substitute  its 
Wj  and  w2  values  in  Equations.  (4.8)  and  (4.9).  The  resulting  values  for  w  and  0 
may  not  correspond  to  any  for  which  Q(w,0)  is  known.  By  virtue  of  equation 
(4.3),  Q(w,0)  will  only  be  known  over  a  uniformly  sampled  set  of  values  for  w 
and  <j>.  In  order  to  determine  Q  at  the  calculated  w  and  <f>,  the  following 
procedure  is  used.  Given  Nw  x  uniformly  located  samples,  Q(w;,<^), 
calculate  a  bilinearly  interpolated  value  of  this  function  at  the  desired  u>  and  <f> 
by  using 


where 


Nw  N, 
i=i  j=i 
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^l^)  ~  |w|<Aw  otherwise, 


(4.10) 


(4.11) 
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A  <f> 

^2(^1  ~  ||  0|  <A0  otherwise;  (4  *2) 

A<f>  and  Au>  are  the  sampling  intervals  for  <j>  and  w,  respectively.  When 
expressed  in  the  manner  shown  above,  bilinear  interpolation  may  be 
interpreted  as  the  output  of  a  filter  whose  impulse  response  is  hjh2. 

The  results  obtained  with  bilinear  interpolation  can  be  considerably 
improved  if  the  sampling  density  in  the  (w,<^)-plane  is  increased  by  using  the 
computationally  efficient  method  of  zero-extending  the  inverse  two-dimensional 
inverse  Fast  Fourier  Transform  (FFT)  of  the  Q(cj;,0j)  matrix.  The  technique 
consists  of  first  taking  a  two-dimensional  inverse  FFT  of  the  Nw  x  matrix 
consisting  of  the  Q(u>j,0j)  values,  zero-extending  the  resulting  Nw  x  array  of 
numbers  to,  perhaps,  mNu  x  nM^  and  then  taking  the  FFT  of  this  new  array. 
The  result  is  an  mn-fold  increase  in  the  density  of  samples  in  the  (u;,0)-plane. 
After  computing  Q(w,^)  at  each  point  of  a  rectangular  grid  by  the  procedure 
outlined  above,  the  object  f(x,y)  is  obtained  by  a  simple  2-D  inverse  FFT. 

The  use  of  bilinear  interpolation  and  zero  padding  are  both  good 
techniques  for  resampling  a  function  but  they  are  used  here  in  a  non  standard 
way.  Typically  interpolation  algorithms  are  derived  assuming  that  the  sampled 
data  can  be  described  as  nearly  linear  (when  using  bilinear  interpolation)  and 
frequency  limited  (when  using  Fourier  domain  zero  padding)  [Con80,  Sto80, 
Act70j.  In  this  application,  when  resampling  the  data  from  a  circular  grid  to  a 
rectangular  grid,  the  function  is  assumed  to  be  smooth  in  the  Fourier  domain. 
This  assumption  is  reasonable  since  the  data  is  assumed  to  be  well  behaved. 

The  interpolation  described  above,  however,  is  carried  out  in  a  rectilinear 
version  of  the  (w,0)  coordinate  system.  Thus  four  points  in  the  ( u,<j) )  space, 
where  data  is  available,  are  first  assumed  to  be  at  the  four  corners  of  a 
rectangle  and  then  the  interpolation  is  calculated  for  a  point  in  the  middle. 
This  is  an  approximation  because  the  four  data  points  actually  define  a  smooth 
function  that  is  defined  along  four  points  on  two  of  the  circular  arcs.  As  will 
be  seen  in  the  reconstructions  the  effect  of  this  approximation  is  small  but  it 
should  be  remembered  when  comparing  interpolation  schemes. 


4.3  Backpropagation  Algorithms 

It  has  recently  been  shown  by  Devaney  [Dev82]  and  Kaveh  [Kav82]  that 
there  is  an  alternative  method  for  reconstructing  images  from  the  diffracted 
projection  data.  This  procedure,  called  the  fdtered-backpropagation  method,  is 
similar  in  spirit  to  the  filtered-backprojection  techniques  which  (due  to  their 
superior  numerical  accuracy)  have  been  one  factor  in  the  enormous  success  of 
x-ray  tomography.  Unfortunately,  whereas  the  filtered-backprojection 
algorithms  also  possess  efficient  implementations,  the  same  cannot  be  said  for 
the  filtered-backpropagation  algorithms.  The  latter  class  of  algorithms  is 
computationally  intensive,  much  more  so  than  the  interpolation  procedure 
discussed  above.  With  regard  to  accuracy,  they  do  not  seem  to  possess  any 
advantage  especially  if  interpolation  is  carried  out  after  increasing  the  sampling 
density  by  appropriate  zero-padding  as  discussed  above. 

The  derivation  of  the  backpropagation  algorithm  will  follow  as  presented 
by  Devaney  [Dev82].  First  consider  the  inverse  Fourier  transform  of  the  object 
function, 


o(rv)  =  777  /  /  rvdK- 

(Z7T)  -00-00 


(4.13) 


This  integral  represents  the  object  function  in  terms  of  the  Fourier  transform 
of  the  object  along  a  rectangular  grid.  As  already  discussed,  a  diffraction 
tomography  experiment  measures  the  Fourier  transform  of  the  object  along 
circular  arcs;  thus  it  will  be  easier  to  perform  the  integration  if  it  is  modified  to 
use  the  projection  data  more  naturally.  This  will  be  done  using  two  coordinate 
transformations:  the  first  one  will  exchange  the  rectangular  grid  for  a  set  of 
semicircular  arcs  and  the  second  will  map  the  arcs  into  their  plane  wave 
decomposition.  ' 

First  exchange  the  rectangular  grid  for  semi-circular  arcs.  To  do  this 
represent  R=(kx,ky)  in  equation  (4.13)  by  the  vector  sum 

R  =  koOHg  (4.14) 


where  s'0=(cos^0,sin^0)  and  T=(cosY,sinX)  are  unit  vectors  representing  the 
direction  of  the  wavevector  for  the  transmitted  and  the  received  plane  waves 
respectively.  This  coordinate  transformation  is  illustrated  in  figure  4.4. 

To  find  the  Jacobian  of  this  transformation  write 


=  k0(sinX-sin^0) 


(4.16) 
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=  k02sin(X-</>0)|  dXd^Q 


=  koy'l-cos^X-^dXd^o 

=  koyi-(s%)2dXd<£0 
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and  then  equation  (4.13)  becomes 
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The  factor  of  —  is  necessary  because  as  discussed  in  section  4.2  the  (X,<^0) 
2 

coordinate  system  gives  a  double  coverage  of  the  (k„,ky)  space. 

This  integral  gives  an  expression  for  the  scattered  field  as  a  function  of  the 
(X,<f>0)  coordinate  system.  While  the  data  that  is  collected  will  actually  be  a 
function  of  <f>0,  the  projection  angle,  and  k,  the  one  dimensional  frequency  of 
the  scattered  field  along  the  receiver  line.  To  make  the  final  coordinate 
transformation  take  the  angle  X  to  be  relative  to  the  (/c,^)  coordinate  system 
diagramed  in  Figure  4.4.  This  is  a  more  natural  representation  since  the  data 
available  in  a  diffraction  tomography  experiment  lies  on  a  semicircle  and 
therefore  the  data  is  available  only  for  0<X<7r.  The  X  integral  in  equation 
(4.20)  above  can  be  rewritten  by  noting 

cosX  =  */k0  (4.21) 


and  therefore 


The  X  integral  becomes 
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Using  the  Fourier  Diffraction  Theorem  as  represented  by  equation  (4.2)  the 
Fourier  transform  of  the  object  function,  O,  can  be  approximated  by  a  simple 
function  of  the  first  order  Born  field,  uR,  at  the  receiver  line.  Thus  the  object 


function  in  equation  (4.24)  can  be  written 

O [k0(sf-^b)]  =  _2'TjUB(/c,'rk0)e  n'°. 


(4.25) 

In  addition  if  a  rotated  coordinate  system  is  used  for  r  =  (£,ff)  where 

£  =  xsin<^  -  ycos <f>  (4.26) 

and 

tj  =  xcos^  +  ysin<£  (4.27) 

then  the  dot  product  ko^-sjj)  can  be  written 

*£  +  (r-  k0)»?-  (4.28) 

The  coordinates  (£,t))  are  illustrated  in  Figure  4.5.  Using  the  results  above  the 

X  integral  is  now  written  as 

ko 

fj-  /  dft|  k\ UB(/£,r-ko)e"hVe+(^  (4.29) 


(4.29) 


and  the  equation  for  the  object  function  in  equation  (4.20)  becomes 

°M  "  7~T7/^^o  /  «j  ^B(K>'pk0)e  J1floeJ'c^+j("rkoK 
(2?r)  o  -k0 


(4.30) 


To  bring  out  the  filtered-backpropagation  implementation,  the  inner 
integration  is  written  here  separately: 


OO 

ry£,>?)  =  —  /  r^(w)H(w)G^(cj)exp(j;<;^)  dw 
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where 
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(4.32) 


r,M  -  UB(K,rk0)e  ,1<0  (4.33) 

Without  the  extra  filter  function  G^w),  the  rest  of  Equation  (4.31)  would 
correspond  to  the  filtering  operation  of  the  projection  data  in  x-ray 
tomography.  The  filtering  as  called  for  by  the  transfer  function  G^(u;)  is  depth 


dependent  due  to  the  parameter  r},  which  is  equal  to  xcos<£  +  ysin^. 

In  terms  of  the  filtered  projections  FI^( ^,17)  in  Equation  (4.31),  the 
reconstruction  integral  of  Equation  (4.30)  may  be  expressed  as 

2ir 

f(x,y)  =  —  J  d^  nJxsin<£  -  ycos <f>  ,  xcos <f>  +  ysin<£)  (4.34) 
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The  computational  procedure  for  reconstructing  an  image  on  the  basis  of 
Equations  (4.31)  and  (4.34)  may  be  presented  in  the  form  of  the  following  steps: 

STEP  1:  In  accordance  with  Equation  (4.31),  filter  each  projection  with  a 
separate  filter  for  each  depth  in  the  image  frame.  For  example,  if 
only  9  depths  are  used  as  shown  in  Figure  4.5,  9  different  filters 
would  need  to  be  applied  to  the  diffracted  projection  shown  there. 
(In  most  cases  for  128  x  128  reconstructive,  the  number  of  discrete 
depths  chosen  for  filtering  the  projection  will  also  be  around  128.  If 
there  are  much  less  than  128,  spatial  resolution  will  be  lost.]  [Cra79j 

STEP  2:  To  each  pixel  (x,y)  in  the  image  frame,  in  accordance  with  Equation 
(4.34)  allocate  a  value  of  the  filtered  projection  that  corresponds  to 
the  nearest  depth  line. 

STEP  3:  Repeat  the  preceding  2  steps  for  all  projections.  As  a  new  projection 
is  taken  up,  add  its  contribution  to  the  current  sum  at  pixel  (x,y). 

The  depth  dependent  filtering  in  Step  1  makes  this  algorithm 
computationally  very  demanding.  For  example,  if  depth  values  are  used, 
the  processing  of  each  projection  will  take  (N^+l)  Fast  Fourier  Transforms 
(FFT’s).  If  the  total  number  of  projections  is  N^,  this  translates  into 
(Nn  +  1)N^  FFT’s.  For  most  N  x  N  reconstructions,  both  and  will  be 
approximately  equal  to  N.  Therefore,  the  filtered-backpropagation  algorithm 
will  require  approximately  N2  FFT’s  compared  to  4N  FFT’s  for  bilinear 
interpolation.  [For  precise  comparisons,  it  must  be  mentioned  that  the  FFT’s 
for  the  case  of  bilinear  interpolation  are  longer  due  to  zero-padding  ] 

Devaney  [I)ev82]  has  also  proposed  a  modified  filtered-backpropagation 
algorithm,  in  which  (’.^(u;)  is  simply  replaced  by  a  single  G^Ju*)  where 
rj0  -  xocos0  +  yo-sind*,  (x0,y0)  being  the  coordinates  of  the  point  where  local 
accuracy  in  reconstruction  is  desired.  (Elimination  of  depth  dependent  filtering 
reduces  the  number  of  FFT’s  to  2 N ^ . ] 


4.4  Signal  Processing  Concerns 

The  reconstruction  algorithms  described  above  and  in  [Pan83]  involve  a 
number  of  signal  processing  steps.  The  following  work  describes  the  quality  of 
the  final  reconstruction  when  small  changes  are  made  to  the  signal  processing 
procedure.  These  changes  are  valid  for  reconstruction  algorithms  using  either 
space  (backpropagation)  or  frequency  domain  interpolation. 

Assuming  first  order  approximations  are  valid,  the  algorithm  for 
reconstructing  an  object  from  diffracted  projections  is  briefly  as  follows: 

1)  Collect  the  data 

2)  Fourier  Transform  each  projection 

3)  Estimate  the  2-dimensional  Fourier  transform  of  the  object  from  the 
transformed  projections 

4)  Perform  a  2-dimensional  inverse  Fourier  transform  to  get  an  estimate  of 
the  object. 

At  each  step  of  this  procedure  signal  processing  theory  suggests  a  number 
of  procedures  to  improve  the  reconstruction.  These  include 

a)  Zero  padding  the  projection  data  to  reduce  the  effects  of  interperiod 
interference.  This  also  increases  the  resolution  in  the  frequency  domain 
and  should  make  interpolation  easier. 

b)  Applying  a  Hamming  window  to  the  projection  data  to  smooth  out  the 
data  at  the  ends  of  the  receiver. 

c)  Multiplying  the  two  dimensional  Fourier  Transform  of  the  object  by  a 
Low  Pass  Filter  (LPF)  (a  Hamming  window  in  this  case)  to  reduce  the 
effects  of  high  frequency  noise. 

These  new  steps  are  illustrated  in  Figure  4.6  where  optional  steps  have  been 
indicated  with  dashed  boxes,  we  have  indicated  the  optional  steps  with  dashed 
boxes. 

To  evaluate  the  effects  of  each  of  these  changes  Figure  4.7  and  Figure  4.8 
shows  the  center  line  of  reconstructions  using  all  eight  possible  combinations  of 
options.  The  data  was  generated  for  a  cylinder  of  radius  IX  and  a  refractive 
index  change  of  .5*^. 

An  important  part  of  the  reconstruction  process  is  filtering  the  projection 
data.  For  efficiency  reasons  the  filter  is  implemented  with  an  FFT  algorithm 
but  these  algorithms  do  not  perform  an  aperiodic  convolution  like  that  used  in 
linear  systems  theory  Instead  a  filter  implemented  with  FFT’s  performs 
circular  convolution. 
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Figure  4.6  The  signal  processing  steps  in  a  typical  diffraction 

tomography  algorithm  are  shown  here.  The  steps  that  are 
needed  are  shown  with  a  solid  box  while  the  optional  steps 
are  shown  with  dashed  lines. 
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Figure  4.8 


The  center  line  of  reconstructions  are  shown  here  with  the 
size  of  the  projection  doubled  and  the  Hamming  window 
added.  All  reconstruction  shown  here  are  with  low  pass 


73 


If  the  data  is  first  padded  with  zeros  so  that  the  new  data  sequence  is 
twice  as  long  as  the  original  then  the  results  produced  by  circular  and  aperiodic 
convolution  are  the  same.  In  addition,  zero  padding  the  original  projection 
increases  the  resolution  of  the  data  in  the  frequency  domain  and  thus  makes  a 
simple  interpolation  scheme  more  accurate.  Unfortunately  the  extra  data  more 
than  doubles  the  computational  time.  (For  example  an  FFT  takes  NIogN 
operations  so  when  N  is  doubled  the  computational  expense  goes  up  by  a  factor 
of 

(2N)log2(2N) 


y&AA 

i _ ( 


Nlog2N 


=  2 


1  + 


log2N 


(4.35) 


Based  on  the  reconstructions  shown  in  Figure  4.7  it  is  possible  to  conclude  that 
doubling  the  size  of  the  projection  data  only  makes  a  small  improvement  in  the 
quality  of  the  reconstructions.  Since  the  extra  zeros  more  than  double  the 
computational  expense  involved  in  filtering  the  data  it  is  probably  best  not  to 
zero  pad. 

A  second  signal  processing  concern  is  due  to  data  truncation.  In  a  real 
world  experiment  it  is  only  possible  to  collect  and  process  a  finite  amount  of 
data.  Generally  this  isn’t  a  problem  since  the  data  eventually  goes  to  zero 
outside  of  some  range  and  the  data  can  be  truncated  without  loss  of 
information.  This  is  certainly  true,  for  example,  in  x-ray  CT  projections  but  is 
not  true  with  diffracted  projections.  With  fields  the  amplitude  decays 

proportional  to  and  consequently  the  projection  data  never  goes  to  zero. 

R 

Mathematically  the  data  truncation  error  can  be  modeled  as  a 
multiplication  in  the  space  domain  by  a  rectangular  window  [Opp75].  In  the 
frequency  domain  this  is  equivalent  to  convolving  the  data  with  a  sine  function 
and  thus  smoot  hes  the  frequency  domain  signal.  A  number  of  windows  like  the 
Hamming  window  have  been  designed  to  reduce  the  effects  of  data  truncation. 

Figure  4.7  shows  that  a  Hamming  window  does  not  have  the  same  positive 
effect  with  diffracted  projections.  In  this  case  most  of  the  high  frequency 
information  is  at  either  end  of  the  projection  and  thus  the  window  only  serves 
to  attenuate  the  high  frequency  components.  This  is  shown  in  Figure  4.9 
where  the  Fourier  transform  of  the  diffracted  projection  is  shown  before  (top 
graph)  and  after  (bottom  graph)  applying  a  Hamming  window  to  the  projection 
data.  The  loss  of  high  frequency  data  caused  by  the  Hamming  window  leads  to 
the  rounded  edges  shown  in  the  reconstructions  shown  in  Figure  4.7. 
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Finally  a  very  large  improvement  is  observed  by  adding  a  Low  Pass  Filter 
before  inverse  Fourier  transforming  the  data.  The  filtering  done  before 
interpolation  includes  a  1/k  term  which  also  serves  to  enhance  the  high 
frequency  noise.  By  adding  a  LPF  this  effect  is  counteracted. 

The  best  results  would  be  obtained  if  it  would  be  possible  to  characterize 
the  spectral  density  of  the  signal  and  the  noise  and  then  design  a  Wiener  filter. 
This  is  a  difficult  problem  and  adequate  results  were  obtained  by  using  a  Low 
Pass  Filter  of  the  form 

m 

w(u>)  =  0.54  +  0.46cos(  — )  (4.36) 

2 

Based  on  the  results  shown  in  Figure  4.7  and  4.8  the  best  reconstructions 
are  obtained  if  a  low  pass  filter  is  used  to  smooth  the  final  reconstruction  but 
that  zero  padding  the  projection  data  and  applying  a  Hamming  window  to  the 
projection  data  do  not  improve  the  results. 

Finally  a  small  improvement  was  made  to  the  backpropagation  algorithm 
by  using  bilinear  interpolation  instead  of  nearest  neighbor.  The 
backpropagation  algorithm  consists  of  both  a  depth  dependent  filter  and  then 


the  addition  to  each  pixel  of  a  portion  of  the  backpropagated  field.  In  the 
original  formulation  each  pixel  is  assigned  the  nearest  neighbor  in  the  field,  but 
as  shown  in  Figure  4.10  even  better  results  are  obtained  if  the  valued  added  to 
each  pixel  is  calculated  using  bilinear  interpolation.  Compared  to  the  expense 
of  doing  the  backpropagation  filter  the  bilinear  interpolation  cost  is 
inconsequential  and  thus  worth  the  effort. 

This  is  also  shown  when  the  Mean  Squared  Error  in  the  reconstructions  is 
computed.  The  table  below  compares  the  error  for  bilinear  ;nterpolation  versus 
nearest  neighbor  and  bilinear  backpropagation.  As  can  be  seen  from  Table  4.1, 
bilinear  interpolation  significantly  improves  the  reconstruction. 


Bilinear  Frequency  Domain  Interpolation 
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Figure  4.10  Reconstructions  of  a  cylinder  are  shown  comparing  nearest 

neighbor  and  bilinear  interpolation  for  frequency  domain  and 
space  domain  reconstruction  algorithms. 
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CHAPTER  6 
LIMITATIONS 

5.1  Introduction 

The  quality  of  a  first  order  diffraction  tomographic  reconstruction  is 
limited  by  both  mathematical  approximations  and  experimental  limitations.  In 
the  derivation  of  a  model  for  the  scattered  fields  either  the  Born  or  the  Rytov 
approximation  is  used  to  solve  the  integral  equations  for  the  scattered  field. 
These  approximations  are  a  source  of  error  and  limit  the  types  of  objects  that 
can  be  imaged  with  diffraction  tomography.  The  only  way  to  reduce  this  type 
of  error  is  to  use  a  better  model  or  a  higher  order  approximation  Better 
models  for  the  scattered  field  will  be  discussed  in  Chapter  6. 

The  experimental  limitations,  on  the  other  hand,  are  caused  by  a  shortage 
of  data.  It  is  only  possible  to  collect  a  finite  amount  of  data  about  the 
scattered  field  and  the  experimental  errors  can  be  attributed  to  interpolation 
errors,  aliasing  and  the  finite  aperture.  Up  to  the  limit  in  resolution  caused  by 
evanescent  waves  and  the  limit  in  quality  due  to  the  Born  and  the  Rytov 
approximations  it  is  possible  to  improve  a  reconstruction  by  collecting 
additional  data. 

Computer  simulations  are  presented  in  this  chapter  illustrating  the  errors 
in  first  order  diffraction  tomography.  To  study  the  effects  of  the  Born  and  the 
Rytov  approximations  it  is  necessary  to  calculate  (or  even  measure)  the  exact 
scattered  fields  and  then  use  the  most  accurate  reconstruction  algorithms 
available.  The  experimental  errors  can  be  minimized  by  calculating  a  large 
number  of  data  points  and  using  a  circularly  symmetric  object  to  reduce  the 
errors  due  to  angular  sampling.  If  experimental  errors  are  minimized  then  the 
only  remaining  source  of  errors  are  caused  by  the  approximations  made  in  the 
reconstruction  algorithm.  As  already  mentioned  the  mathematical  limitations 
on  the  reconstructions  are  a  function  of  the  object's  size  and  refractive  index. 

The  experimental  errors  are  highlighted  by  minimizing  the  algorithmic 
errors.  This  can  be  done  in  two  ways.  The  more  straightforward  method  is  to 
choose  a  small  object  with  a  small  change  in  refractive  index.  As  the  size  of 
the  object  and  its  refractive  index  are  reduced  both  the  Born  and  the  Rytov 


approximations  more  accurately  model  the  exact  scattered  field.  A  second 
approach  is  to  assume  the  Born  or  Rytov  approximations  are  valid  then  use  the 
Fourier  Diffraction  Theorem  to  generate  the  scattered  fields  from  the  Fourier 
transform  of  the  object.  In  both  cases  the  amount  of  data  calculated  is  varied 
to  highlight  the  different  experimental  errors. 

5.2  Mathematical  Limitations 

In  diffraction  tomography  there  are  different  approximations  involved  in 
the  forward  and  inverse  directions.  In  the  forward  process  it  is  necessary  to 
assume  that  the  object  is  weakly  scattering  so  that  either  the  Born  or  the 
Rytov  approximations  can  be  used.  Once  an  expression  for  the  scattered  field 
is  derived  it  is  necessary  to  not  only  to  measure  the  scattered  fields  but  then 
numerically  implement  the  inversion  process. 

To  study  the  limits  of  the  mathematical  approximations  the  exact  field  for 
the  scattered  field  from  a  cylinder  as  shown  by  Weeks  [Wee64]  and  by  Morse 
and  Ingard  (Mor6H]  is  calculated  for  cylinders  of  various  sizes  and  refractive 
index  In  the  simulates  that  follow  a  single  plane  wave  of  unit  wavelength  is 
incident  on  the  cylinder  and  the  scattered  field  is  measured  along  a  line  at  a 
distance  of  100  wavelengths  from  the  origin.  In  addition  all  refractive  index 
changes  are  modeled  as  monopole  (omnidirectional)  scatterers.  By  doing  this 
the  directional  dependence  of  dipole  scatterers  does  not  have  to  be  taken  into 
account . 

At  the  receiver  line  the  received  wave  is  measured  at  512  points  spaced  at 
1/2  wavelength  intervals.  In  all  cases  the  rotational  symmetry  of  a  single 
cylinder  at  the  origin  is  used  to  reduce  the  computation  time  of  the 
simulations.  Since  all  projections  are  identical  this  eliminates  any  angular 
interpolation  error. 


5.2.1  Evaluation  of  the  Born  Approximation 


In  using  the  Born  approximation  it  is  necessary  to  assume  that  the 
amplitude  of  the  scattered  field  is  small  compared  to  the  incident  field.  As 
already  discussed  the  Born  approximation  is  most  sensitive  to  phase  changes 
and  this  will  be  shown,  first  qualitatively  and  then  quantitatively.  From 
Chapter  2  the  phase  change  is  given  by 


Phase  Change  =  47m 


6  X 


(5.1) 


where  the  cylinder  has  a  radius  of  ‘a’  and  a  refractive  index  of  1  +n5. 


The  results  shown  in  Figure  5.1  are  for  cylinders  of  four  different  refractive 
indices.  In  addition  Figure  5.2  shows  plots  of  each  reconstruction  along  a  line 
through  the  center  of  the  cylinder.  Notice  that  the  y  coordinate  (refractive 
index)  of  the  center  line  is  plotted  in  terms  of  change  from  unity. 

Simulations  are  shown  for  refractive  indices  that  range  from  .1%  change 
(refractive  index  of  1.001)  to  a  40%  change  (refractive  index  of  1.4).  For  each 
refractive  index,  cylinders  of  size  1,  2,  4  and  10  wavelengths  are  shown.  This 
gives  a  range  of  phase  changes  across  the  cylinder  (see  equation  (5.1)  above) 
from  .004jt  to  16jt. 

Clearly,  all  the  cylinders  of  refractive  index  1.001  in  Figure  5.1  are 
perfectly  reconstructed.  As  equation  (5.1)  predicts  the  results  get  worse  as  the 
product  of  refractive  index  and  radius  gets  larger.  The  largest  refractive  index 
that  is  successfully  reconstructed  is  for  the  cylinder  in  Figure  5.1  of  radius  1 
wavelength  and  a  refractive  index  that  differed  by  20%  from  the  surrounding 
medium. 

While  it  is  difficult  to  evaluate  quantitatively  the  three  dimensional  plots 
it  is  certainly  reasonable  to  conclude  that  only  cylinders  where  the  phase 
change  across  the  object  is  less  than  or  equal  to  .Sn  are  adequately 
reconstructed.  In  general  the  reconstruction  for  each  cylinder  where  the  phase 
change  across  the  cylinder  is  greater  than  ir  shows  severe  artifacts  near  the 
center.  This  limitation  in  the  phase  change  across  the  cylinder  is  consistent 
with  the  condition  described  in  Chapter  2. 

5.2.2  Evaluation  of  the  Rytov  Approximation 

Figure  5.3  shows  the  simulated  results  for  16  reconstructions  using  the 
Rytov  approximation.  To  emphasize  the  insensitivity  of  the  Rytov 
approximation  to  large  objects  the  largest  object  simulated  has  a  diameter  of 
100X. 

It  should  be  pointed  out  that  the  rounded  edges  of  the  IX  reconstructions 
are  not  due  to  any  limitation  of  the  Rytov  approximation  but  instead  are  the 
result  of  a  two  dimensional  low  pass  filtering  of  the  reconstructions.  Recall 
that  for  a  transmission  experiment  an  estimate  of  the  object’s  Fourier 
transform  is  only  available  up  to  frequencies  less  than  v/2k0.  Thus  the 
reconstructions  shown  in  Figure  5.3  show  the  limitations  of  both  the  Rytov 
approximation  and  the  Fourier  Diffraction  Theorem. 
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Figure  5.3 


Simulated  reconstructions  using  the  Born  approximation  for 
16  objects  with  four  refractive  indices  between  1.001  and  1.10 
and  four  radii  between  1  and  100X. 


V*  t  *  «  *  «  '  i  * 

>»  *  '  .  fj 


•  .  •  ,•/«>>  .  •  /  V  -  •  *.•  V/vV  .•  V.  .•  /.  w*.  . 


88 


t 


5.2.3  Comparison  of  the  Born  and  Rytov  Approximation 

Reconstructions  using  exact  scattered  data  show  the  similarity  of  the  Born 
and  the  Rytov  approximation.  Within  the  limits  of  the  Fourier  Diffraction 
Theorem  the  reconstructions  in  Figure  5.1  and  5.3  of  a  IX  object  with  a  small 
refractive  index  are  similar.  In  both  cases  the  reconstructed  change  in 
refractive  index  is  close  to  that  of  the  simulated  object. 

The  two  approximations  differ  for  objects  that  have  a  large  refractive 
index  change  or  have  a  largA  radius.  The  Born  reconstructions  are  good  at  a 
large  refractive  index  as  long  as  the  phase  shift  of  the  incident  field  as 
predicted  by  equation  (5.1)  is  less  than  ir. 

On  the  other  hand  the  Rytov  approximation  is  very  sensitive  to  the 
refractive  index  but  produce  excellent  reconstructions  for  objects  as  large  as 
100X.  Unfortunately  for  object  with  a  refractive  index  larger  than  \%  the 
Rytov  approximation  quickly  deteriorates. 

In  addition  to  the  qualitative  studies  a  quantitative  study  of  the  error  in 
the  Born  and  Rytov  reconstructions  is  also  possible.  As  a  measure  of  error  the 
relative  mean  squared  error  in  the  reconstruction  of  the  object  function  is 
integrated  over  the  entire  plane.  If  the  actual  object  function  is  o(r)  and  the 
reconstructed  object  function  is  o' (r)  then  the  relative  Mean  Squared  Error 
(MSE)  is  given  by 

JJ[o(r)-o'(r)]2dr  ,  4 

//!°(r)]2  •  '  •*’ 

This  study  presents  the  Mean  Squared  Error  for  120  reconstructions  based 
on  the  exact  scattered  fields  from  a  cylinder.  In  each  case  a  512  point  receiver 
line  is  at  a  distance  of  10  X  from  the  center  of  the  cylinder.  Both  the  receiver 
line  and  the  object  reconstruction  are  sampled  at  1/4  X  intervals. 

The  plots  of  Figure  5.4  present  a  summary  of  the  mean  squared  error  for 
cylinders  of  1,  2  and  3  X  in  radius  and  for  twenty  refractive  indices  between 
1.01  and  1.20.  In  each  case  the  error  for  the  Born  approximation  is  shown  as  a 
solid  line  while  the  Rytov  reconstruction  is  shown  as  a  dashed  line.  The  data 
used  for  these  simulations  was  the  exact  scattered  fields  from  a  cylinder 
measured  at  512  receiver  points  along  a  receiver  line  10X  from  the  center  of  the 
cylinder. 

Many  researchers  [Kav82,  Kel69,  Sou83j  have  postulated  that  the  Rytov 
approximation  is  clearly  superior  to  the  Bom  but  as  the  actual  reconstructions 
represented  by  Figure  5.4  show  for  a  IX  cylinder  this  is  not  necessarily  true.  It 
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The  relative  mean  squared  error  for  reconstructions  using  the 
Born  (solid  line)  and  Rytov  (dashed  line)  approximations. 
The  error  for  a  total  of  60  objects  with  a  radius  of  1,  2  and  3 
wavelengths  are  shown. 


is  interesting  to  note  that,  while  the  Rytov  approximation  shows  a  steadily 
increasing  error  with  higher  refractive  indices,  the  error  in  the  Born 
reconstruction  is  relatively  constant  until  a  threshold  is  reached.  For  the  2X 
and  3X  cylinders,  this  breakpoint  occurs  at  a  phase  shift  across  the  cylinder  of 
0.6  and  0.7jt.  Thus,  a  criteria  for  the  validity  of  the  Born  approximation  is 
that  the  product  of  the  radius  of  the  cylinder  in  wavelengths  and  the  change  in 
refractive  index  must  be  less  than  .175. 

Figure  5.5  presents  a  summary  of  the  relative  mean  squared  errors  for 
cylinders  with  refractive  indices  of  1.01,  1.02  and  1.03  for  forty  radii  between  1 
and  40X.  Because  the  size  of  the  cylinders  varied  by  a  factor  of  forty,  the 
simulations  parameters  were  adjusted  accordingly.  For  a  cylinder  of  radius  R, 
the  scattered  field  was  calculated  for  512  receivers  along  a  line  2R  from  the 
center  of  the  cylinder  and  spaced  at  1/16R  intervals. 

In  each  of  the  simulations,  the  Born  approximation  is  only  slightly  better 
than  the  Rytov  approximation  until  the  Born  approximation  crosses  its 
threshold  with  a  phase  shift  of  0.7jt.  Because  the  error  in  the  Rytov 
approximation  is  relatively  flat,  it  is  clearly  superior  for  large  objects  with 
small  refractive  indices.  Using  simulated  data  and  the  Rytov  approximation 
objects  as  large  as  2000X  in  radius  have  been  reconstructed. 

5.3  Evaluation  of  Reconstruction  Algorithms 

To  study  the  approximations  involved  in  the  reconstruction  process  it  is 
necessary  to  calculate  scattered  data  assuming  the  forward  approximations  are 
valid.  This  can  be  done  in  one  of  two  different  ways.  As  already  discussed  the 
Born  and  Rytov  approximations  are  valid  for  small  objects  and  small  changes 
in  refractive  index.  Thus,  if  the  exact  scattered  field  for  a  small  and  weakly 
scattering  object  is  calculated  then  it  can  be  assumed  that  either  the  Born  or 
the  Rytov  approximations  is  exact. 

A  better  approach  is  to  recall  the  Fourier  Diffraction  Theorem,  which  says 
that  the  Fourier  transform  of  the  scattered  field  is  proportional  to  the  Fourier 
transform  of  the  object  along  a  circular  arc.  Since  this  theorem  is  the  basis  for 
the  first  order  inversion  algorithms  if  it  is  assumed  correct  the  approximations 
involved  in  the  reconstruction  process  can  be  studied. 

If  the  Fourier  Diffraction  Theorem  holds,  the  exact  scattered  field  can  be 
calculated  exactly  for  objects  that  can  be  modeled  as  ellipses.  The  analytic 
expression  for  the  Fourier  transform  of  the  object  along  an  arc  is  then 
proportional  to  the  scattered  fields.  This  procedure  is  fast  and  allows  the 
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Figure  5.5  The  relative  mean  squared  error  for  reconstructions  using  the 

Born  (solid  line)  and  Rytov  (dashed  line)  approximations. 
The  error  for  a  total  of  120  objects  with  a  refractive  index  of 
1.01,  1.02  and  1.03  are  shown. 


scattered  fields  to  be  calculated  for  testing  reconstruction  algorithms  and 
experimental  parameters. 

To  illustrate  the  accuracy  of  the  interpolation  based  algorithms,  the  image 
in  Figure  5.6a  will  be  used  as  a  test  “object”  for  showing  computer  simulation 
results.  Figure  5.6a,  with  its  gray  levels  as  shown  in  Figure  5.6b,  is  a 
modification  of  the  Shepp  and  Logan  “phantom”  described  in  [She74]  to  the 
case  of  diffraction  imaging  [Dev83,  Pan83).  The  gray  levels  shown  in  Figure  5.6 
represent  the  refractive  index  values.  This  test  image  is  a  superposition  of 
ellipses,  with  each  ellipse  being  assigned  a  refractive  index  value  as  shown  in 
Table  5.1. 


Table  5.1.  Summary  of  parameters  for  diffraction 
tomography  simulations 


Center 

Coordinate 

Major 

Axis 

Minor 

Axis 

Rotation 

Angle 

Refractive 

Index 

(0,0) 

0.92 

0.69 

90 

1.0 

(0,-0.0184) 

0.874 

0.6624 

90 

-0.5 

(0.22,0) 

0,31 

0.11 

72 

-0.2 

(-0.22,0) 

0.41 

0.16 

108 

f 

© 

(0,0.35) 

0.25 

0.21 

90 

0.1 

(0,0.1) 

0.046 

0.046 

0 

0.15 

(0,-0. 1) 

0.046 

0.046 

0 

0.15 

(-0.08,-0.605) 

0.046 

0.023 

0 

0.15 

(0,-0.605) 

0.023 

0.023 

0 

0.15 

(0.06,-0.605) 

0.046 

0.023 

90 

0.15 

A  major  advantage  of  using  an  image  like  Figure  5.6a  for  computer 
simulation  is  that  the  analytical  expressions  for  the  transforms  of  the  diffracted 
projections  can  be  written.  The  Fourier  transform  of  an  ellipse  of  semi-major 
and  semi-minor  axes  A  and  B,  respectively,  is  given  by 
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Figure  5.6  A  modified  version  of  the  Shepp  and  Logan  head  phantom  is 

used  to  test  reconstruction  algorithms.  The  numbers 
represent  the  relative  change  in  refractive  index  from  the 
background  value  of  1.0. 


where  u  and  v  are  spatial  angular  frequencies  in  the  x-  and  y-directions, 
respectively;  and  Jj  is  a  Bessel  function  of  the  first  kind  and  order  1.  When  the 
center  of  this  ellipse  is  shifted  to  the  point  (xj.y,),  and  the  angle  of  the  major 
axis  tilted  by  a,  as  shown  in  Figure  5.7b,  its  Fourier  transform  becomes 


-j(ux,  +  vyi)  . 


27rAJ  i  B  ((u')A/B)2  +  (v')2 


(K  )A/B)2  +  (v'  )2 


where 


u'  =  u  cos  a  +  v  sin  a 


v'  =  -u  sin  a  +  v  cos  a.  (5.6) 

Now  consider  the  situation  where  the  ellipse  is  illuminated  by  a  plane.  By 
the  Fourier  Diffraction  Theorem  the  Fourier  transform  of  the  transmitted  wave 
fields  measured  on  the  receiver  line  is  given  by  the  values  of  the  above  function 
on  a  circular  arc.  For  the  test  object  of  Figure  5.7b,  if  weak  scattering  is 
assumed  and  therefore  there  is  no  interaction  among  the  ellipses,  the  Fourier 
transform  of  the  total  forward  scattered  field  measured  on  a  line,  is  the  sum  of 
the  values  of  functions  like  (5.4)  over  the  circular  arc.  This  procedure  is  used 
to  generate  the  diffracted  projection  data  for  the  test  image. 

It  should  be  mentioned  that  by  using  this  procedure  to  calculate  the 
diffracted  projection  data  only  the  accuracy  of  the  reconstruction  algorithm  is 
tested,  without  checking  whether  or  not  the  “test  object ’’  satisfies  the  underlying 
assumption  of  weak  scattering.  In  order  to  test  this  crucial  assumption,  it  is 
necessary  to  generate  the  forward  scattered  data  of  the  object.  For  multi- 
component  objects,  such  as  the  one  shown  in  Figure  5.6a,  it  is  very  difficult  to 
do  so  due  to  the  interactions  between  the  components. 

Pan  and  Kak  [Pan83]  presented  the  simulations  shown  in  Figure  5.8. 
Using  a  combination  of  increasing  the  sampling  density  by  zero  padding  the 
signal  and  bilinear  interpolation,  results  are  obtained  in  2  minutes  of  CPU  time 
on  a  VAX  11/780  minicomputer  with  a  floating  point  accelerator  (FPA).  The 
reconstruction  is  shown  over  a  128  by  128  grid  using  64  views  and  128  receiver 
positions.  The  number  of  operations  required  to  carry  out  the  interpolation 
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Figure  5.7  Assuming  the  Fourier  Slice  Theorem  is  valid  the  scattered 

field  can  easily  be  computed  as  the  values  of  the  Fourier 
transform  of  the  rotated  and  translated  ellipse. 


Figure  5.8  The  above  images  show  the  results  of  using  the  (al 

interpolation,  (b)  backpropagation  and  (c)  modifiea 
backpropagation  algorithms  on  reconstruction  quality.  The 
solid  lines  in  the  graphs  represent  the  reconstructed  value 
along  a  line  through  the  center  of  the  three  ellipses  at  the 
bottom  of  the  phantom. 


and  invert  the  object  function  is  on  the  order  of  N2iogN.  The  resulting 
reconstruction  is  shown  in  Figure  5.8a. 

Figure  5.8b  represents  the  result  of  back  propagating  the  data  to  128 
depths  for  each  view  while  Figure  5.8c  is  the  result  of  back  propagation  to  only 
a  single  depth  centered  near  the  three  small  ellipses  at  the  bottom  of  the 
picture.  The  results  were  calculated  on  a  VAX  11/780  minicomputer  with  a 
Floating  Point  Accelerator  (FPA)  and  the  resulting  reconstructions  were  done 
over  a  128  by  128  grid.  Like  the  previous  image  the  input  data  consists  of  64 
projections  of  128  points  each. 

There  is  a  significant  difference  in  not  only  the  reconstruction  time  but 
also  the  resulting  quality.  While  the  modified  back  propagation  only  took  1.25 
minutes  the  resulting  reconstruction  is  much  poorer  than  the  full  back 
propagation  which  took  30  minutes  of  CPU  time.  A  comparison  of  the  various 
algorithms  is  shown  in  Table  5.2. 


Table  5.2.  Comparison  of  Algorithms 


Algorithm 

Complexity 

CPU  Time 

Interpolation 

N2logN 

2  Minutes 

Back  Propagation 

NdN*NlogN 

30  Minutes 

Modified  Back  Propagation 

N,NlogN 

1.25  Minutes 

6.4  Experimental  Limitations 

In  addition  to  the  limits  on  the  reconstructions  imposed  by  the  Born  and 
the  Rytov  approximations  there  are  also  experimental  limitations.  These 
additional  factors  are  caused  by 

•  Wave  propagation  in  free  space 

•  Sampling  the  data  along  the  receiver  line 

•  Finite  receiver  length 

•  Limited  views  of  the  object 

In  inverse  scattering  theory  the  measured  complex  amplitude  of  a  received 
wave  is  sampled,  filtered  and  then  interpolated  to  estimate  the  Fourier 
transform  of  the  object  function.  The  reconstruction  process  is  linear  because 


it  consists  only  of  filtering  the  data  and  then  calculating  its  inverse  Fourier 
transform. 

The  first  three  factors  each  can  be  modeled  as  a  simple  constant  low-pass 
filtering  of  the  scattered  field.  Because  the  entire  process  is  linear  the  net  effect 
is  a  single  low  filter  at  the  lowest  of  the  three  frequencies.  The  experiment  can 
be  optimized  by  adjusting  the  parameters  of  the  experiment  so  that  each  low 
pass  filter  cuts  off  at  the  same  frequency. 

A  limited  number  of  views  also  can  be  modeled  as  a  low  pass  filter.  In  this 
case,  though,  the  cutoff  frequency  varies  with  the  radial  direction. 

6.4.1  Evanescent  Waves 

The  most  fundamental  limitation  is  that  evanescent  waves  are  ignored. 
Since  these  waves  have  a  complex  wavenumber  they  are  severely  attenuated 
over  a  distance  as  small  as  several  wavelengths.  This  limits  the  highest 
received  wavenumber  to 

k„„  =  y  (5.7) 

This  is  a  fundamental  limit  of  the  propagation  process  and  the  resolution  of  the 
reconstruction  can  only  be  improved  by  moving  the  experiment  to  a  higher 
frequency  (or  shorter  wavelength.) 

6.4.2  Sampling  the  Received  Wave 

After  the  wave  has  been  scattered  by  the  object  and  propagated  to  the 
receiver  line  it  must  be  measured.  This  is  usually  done  with  the  a  point 
receiver.  Unfortunately  it  is  not  possible  to  sample  at  every  point,  so  a  non 
zero  sampling  interval  must  be  chosen.  This  introduces  a  measurement  error 
into  the  process.  By  the  Nyquist  theorem  this  can  be  modeled  as  a  low  pass 
filtering  operation,  where  the  highest  measured  frequency  is  given  by 

km.„  =  J  (5.8) 

and  T  is  the  sampling  interval.  Of  course  this  analysis  has  ignored  the  non¬ 
linear  effects  of  aliasing  and  the  resulting  frequency  shifts  that  occur. 


5.4.3  The  Effects  of  a  Finite  Receiver  Length 

Not  only  are  there  physical  limitations  on  the  finest  sampling  interval  but 
usually  there  is  a  limitation  on  the  amount  of  data  that  can  be  collected.  This 
generally  means  that  samples  of  the  received  waveform  will  be  collected  at  only 
a  finite  number  of  points  along  the  receiver  line.  This  is  usually  justified  by 
taking  data  along  a  line  long  enough  so  that  the  unmeasured  data  can  be  safely 
ignored.  Because  of  the  wave  propagation  process  this  also  introduces  a  low 
pass  filtering  of  the  received  data. 

Consider  for  a  moment  a  single  scatterer  at  some  distance,  10,  from  the 
receiver  line.  The  wave  propagating  from  this  single  scatterer  is  a  cylindrical 
wave  in  two  dimensions  or  a  spherical  wave  in  three  dimensions.  This  effect  is 
diagramed  in  Figure  5.9.  It  is  easy  to  see  that  the  spatial  frequencies  vary  with 
the  position  along  the  receiver  line. 

An  optician  studying  this  problem  would  be  interested  in  knowing  the 
resolving  power  of  the  system  as  a  function  of  the  size  of  the  aperture  [G0068]. 
The  analysis  would  normally  be  carried  out  assuming  that  the  object  is  far 
enough  from  the  aperture  so  that  it  can  be  assumed  it  is  in  the  aperture’s  far 
field.  But  in  this  work  the  frequency  content  of  the  measured  field  is  a  limiting 
factor  in  the  reconstruction  quality  so  the  effect  of  a  limited  aperture  will  be 
analyzed  with  an  emphasis  on  the  spatial  frequency  content  of  the  received 
field.  The  two  approaches  to  be  considered  here  use  a  point  source  and  analyze 
the  frequency  content  at  the  aperture.  Since  all  points  in  space  are  in  the  far 
field  of  a  point  source  this  analysis  gives  identical  results  to  classical  optics 
theory. 

It  is  easier  to  analyze  the  effect  by  considering  the  expanding  wave  to  be 
locally  planar  at  any  point  distant  from  the  scatterer.  At  the  point  on  the 
receiver  line  closest  to  the  scatterer  there  is  no  spatial  variation  [Goo68j.  This 
corresponds  to  receiving  a  plane  wave  or  a  received  spatial  frequency  of  zero. 

Higher  spatial  frequencies  are  received  at  points  along  the  receiver  line 
that  are  farther  from  the  origin.  The  received  frequency  is  a  function  of  the 
sine  of  the  angle  between  the  direction  of  propagation  and  a  perpendicular  to 
the  receiver  line.  This  function  is  given  by 

My)  =  kmMsin0  (5-°) 

where  9  is  the  angle  and  kmM  is  the  wavenumber  of  the  incident  wave.  Thus 
at  the  origin,  the  angle,  9,  is  zero  and  the  received  frequency  is  zero.  Only  at 
infinity  does  the  angle  become  equal  to  ninety  degrees  and  the  received  spatial 
frequency  approach  the  theoretical  maximum. 


This  reasoning  can  be  justified  on  a  more  theoretical  basis  by  considering 
the  phase  of  the  propagating  wave.  The  received  Wave  at  a  point  (x=l0y)  due 
to  a  scatterer  at  the  origin  is  given  by 


u(x=lo,y) 


ik0\/x2+y4 


(5.10) 


w  '  Vx2+y2 

The  instantaneous  spatial  f'^quency  along  the  receiver  line  (y  varies)  of  this 
wave  can  be  found  by  taking  the  partial  derivative  of  the  phase  with  respect  to 
y  [Gag78] 

phase  =  k0Vy2+x2  (5.11) 


k  =  k°y 
recv 


(5.12) 


where  krecv  is  the  spatial  frequency  received  at  the  point  (x=lo,y).  From  Figure 
5.9  it  is  easy  to  see  that 


sin#  = 


(5.13) 


Vx2+y2 

and  therefore  equation  (5.9)  and  (5.12)  are  equivalent. 

This  relation,  (5.12),  can  be  inverted  to  give  the  length  of  the  receiver  line 
for  a  given  maximum  received  frequency,  kmM.  This  becomes 


y  =  ±- 


k  x 

•'max 


A^-k  2 

*0  ’'max 


(5.14) 


Since  the  highest  received  frequency  is  a  monotonically  increasing  function 
of  the  length  of  the  receiver  line  it  is  easy  to  see  that  by  limiting  the  sampling 
of  the  received  wave  to  a  finite  portion  of  the  entire  line  that  a  low  passed 
version  of  the  entire  scattered  wave  is  measured.  The  highest  measured 
frequency  is  a  simple  function  of  the  distance  of  the  receiver  line  from  the 
scatterer  and  the  length  of  measured  data.  This  limitation  can  be  better 
understood  if  the  maximum  received  frequency  is  written  as  a  function  of  the 
angle  of  view  of  the  receiver  line.  Thus  substituting 


tan(0)  = 


(5.15) 


it  is  easy  to  see  that 
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(5.16) 


*  +  12 

x 


^recv 


(5.17) 


reCV  v/tan^  +  1 

Thus  krecv  is  a  monotonically  increasing  function  of  the  angle  of  view,  6.  It  is 
easy  to  see  that  the  maximum  received  spatial  frequency  can  be  increased  by 
either  moving  the  receiver  line  closer  to  the  object  or  by  increasing  the  length 
of  the  receiver  line. 


6.4.4  Evaluation  of  the  Experimental  Effects 

Exact  scattered  data  is  used  to  verify  the  optimum  experimental  values 
and  the  effect  of  a  finite  receiver  length  is  shown  in  Figure  5.10.  The  spatial 
frequency  content  of  a  wave  is  found  by  taking  the  FFT  of  the  sampled  points 
along  the  receiver  line  and  is  compared  to  the  theoretical  result  as  predicted  by 
the  Fourier  Transform  of  the  object.  The  theory  predicts  that  more  of  the 
high  frequency  components  will  be  present  as  the  length  of  the  receiver  line 
increases  and  this  is  confirmed  by  simulation. 

While  the  above  derivation  only  considered  a  single  scatterer  it  is  also 
approximately  true  for  many  scatterers  collected  at  the  origin.  This  is  possible 
because  the  inverse  reconstruction  process  is  linear  and  each  point  in  the  object 
scatterers  an  independent  cylindrical  wave. 


5.4.5  Optimiiation 

Since  each  of  these  three  factors  is  independent,  their  effect  in  the 
frequency  domain  can  be  found  by  simply  multiplying  each  of  their  frequency 
responses  together.  As  has  been  described  above  each  of  these  effects  can  be 
modeled  as  a  simple  low-pass  filter  of  the  theoretical  data  so  the  combined 
effect  is  also  a  low  pass  filter  but  at  the  lowest  frequency  of  the  cut-off  of  the 
three  effects. 

First  consider  the  effect  of  ignoring  the  evanescent  waves.  Since  the 
maximum  frequency  of  the  received  wave  is  limited  by  the  propagation  filter  to 


*  V"  A1  *.*  *  *  »  *  i  ‘  »  *  *  »  *  •  *  *  "  v'  »  *  »  "  k  '  «  *  i  \  <»  »  •  v  \  *■  v  -  »  • 

'•.•vVA'V-.v :■ .  y -v •*- 


yfr.v 


V  /  / 

Vc.  ‘  i 

v  v  V 
*  ■»  < 

v 

•> 


■ivKM 


;v,v/ 


V  %  ■ 
,■*  V  %• 

’'w  N.  ^ 


>.v 


-V 


*«\*%«  *. 


t'Vy. 


IU 


V  •  . 

- 


it  is  easy  to  combine  this  expression  with  the  expression  for  the  Nyquist 
frequency  into  a  single  expression  for  the  smallest  “interesting”  sampling 
interval.  This  is  given  by 

—  ^meu  (5.19) 


Y  =  t  (520) 

and  rearranging 

T>1/2X.  (5.21) 

Thus  if  the  received  waveform  is  sampled  with  a  sampling  interval  of  more 
than  1/2  wavelength  the  measured  data  might  not  be  a  good  estimate  of  the 
received  waveform  because  of  aliasing.  On  the  other  hand  it  is  not  necessary  to 
sample  the  received  waveform  any  finer  than  1/2  wavelength  since  this 
provides  no  additional  information.  Doing  this  it  is  possible  to  conclude  that 
the  sampling  interval  should  be  greater  than  1/2  wavelength. 

In  general  the  experiment  will  also  be  constrained  by  the  number  of  data 
points  (M)  that  can  be  measured  along  the  receiver  line.  The  distance  from  the 
object  to  the  receiver  line  will  be  considered  a  constant  in  the  derivation  that 
follows.  If  the  received  waveform  is  sampled  uniformly,  the  range  of  the 
receiver  line  is  given  uniquely  by 

y„,„  =  (5.22) 

This  is  also  shown  in  Figure  5.9. 

For  a  receiver  line  at  a  fixed  distance  from  the  object  and  a  fixed  number 
of  receiver  points  this  is  a  classical  optimization  problem.  As  the  sampling 
interval  is  increased  the  length  of  the  receiver  line  increases  and  more  of  the 
received  wave’s  high  frequencies  are  measured.  On  the  other  hand  increasing 
the  sampling  interval  lowers  the  maximum  frequency  that  can  be  measured 
before  aliasing  occurs. 

The  optimum  value  of  T  can  be  found  by  setting  the  cutoff  frequencies  for 
the  Nyquist  frequency  equal  to  the  highest  received  frequency  due  to  the  finite 
receiver  length  and  then  solving  for  the  sampling  interval.  If  this  constraint  is 
not  met  then  some  of  the  information  that  is  passed  by  one  process  will  be 
attenuated  by  the  others.  This  results  in 


evaluated  at 


k0= 


2?r 


and 


y  - 


MT 


Solving  for  T2  the  optimum  value  for  T  is  given  by 
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64 


+  M2  +  M 
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Making  the  substitution 


«  — 


XM 


the  optimum  sampling  interval  is  given  by 

_  \/64q2  +  1  +  1 


X 
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(5.24) 


(5.25) 


(5.26) 


(5.27) 


(5.28) 


This  substitution  is  similar  to  the  tan#  substitution  that  is  made  in  the 
heuristic  approach  above.  Also  notice  that  the  smallest  positive  value  that  the 
sampling  interval  can  become  is  1/2  wavelength.  This  corresponds  to  the 
Nyquist  frequency  for  a  propagating  wave. 

The  optimum  sampling  interval  is  confirmed  by  simulations.  Again,  using 
the  method  described  above  for  calculating  the  exact  scattered  fields,  four 
simulations  are  shown  of  an  object  of  radius  10  wavelengths  using  a  receiver 
that  is  100  wavelengths  from  the  object.  In  each  case  the  number  of  receiver 
positions  is  fixed  at  64.  The  resulting  reconstructions  for  a  sampling  interval  of 
.5,  1,  1.5  and  2  wavelengths  are  shown  in  Figure  5.10.  Equation  (5.28)  predicts 
an  optimum  sampling  interval  of  1.3  wavelengths  and  this  is  confirmed  by  the 
simulations.  The  best  reconstruction  occurs  with  a  sampling  interval  between  1 
and  1.5  wavelengths. 
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6.4.6  Limited  Views 

In  many  applications  it  is  not  possible  to  generate  or  receive  plane  waves 
from  all  directions.  The  effect  of  this  it  to  leave  holes  where  there  is  no 
estimate  of  the  Fourier  transform  of  the  object. 

Since  the  ideal  reconstruction  algorithm  produces  an  estimate  of  the 
Fourier  transform  of  the  object  for  all  frequencies  within  a  disk  a  limited 
number  of  views  introduces  a  selective  filter  for  areas  where  there  is  no  data. 
As  shown  by  Devaney  [Dev84]  for  the  VSP  case  a  limited  number  of  views 
degrades  the  reconstruction  by  low  pass  filtering  the  image  in  certain 
directions.  Devaney’s  results  are  reproduced  in  Figures  5.11  and  5.12. 
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CHAPTER  6 

HIGHER  ORDER  APPROXIMATIONS 
TO  THE  SCATTERED  FIELD 


I 

i 

) 

6.1  Introduction 

Reconstruction  algorithms  based  on  the  theory  presented  in  Chapter  2 
generate  severe  artifacts  for  objects  larger  than  a  few  wavelengths  or  a  change 
in  refractive  index  greater  than  a  few  percent.  These  reconstruction  algorithms 
»  are  limited  by  the  first  Born  or  the  first  Rytov  approximations  and  thus  any 

improvement  in  the  reconstructions  will  be  accomplished  by  modeling  the 
scattered  fields  more  accurately.  With  this  more  accurate  model  it  will  again 
be  necessary  to  invert  the  equation  to  arrive  at  an  estimate  of  the  object. 

This  section  will  describe  iterative  techniques  that  more  accurately  model 
the  scattered  fields.  In  addition  to  the  theoretical  derivations  a  number  of 
simulations  will  be  presented  and  the  convergence  of  each  series  will  be 
discussed. 

Two  approaches  to  more  accurately  compute  the  scattered  field  will  be 
discussed  here.  Both  approaches  are  iterative  but  they  differ  widely  in  their 
philosophy.  The  first  approach,  known  as  a  fixed  point  method,  assumes  an 
initial  guess  that  is  a  small  perturbation  from  the  correct  solution.  If  the 
scattering  integral  is  a  “contracting”  operator  then  this  procedure  converges  to 
the  correct  answer.  The  second  approach,  based  on  a  matrix  formulation  of 
the  discrete  wave  equation,  refines  an  initial  guess  by  projecting  it  onto 
hyperplanes.  Since  each  projection  reduces  the  error  this  method  always 
converges.  Unfortunately  the  matrix  formulation  is  more  expensive  (requires 
more  calculations)  than  the  fixed  point  method  by  a  factor  of  several  hundred. 


6.2  The  Singularity  of  the  Green’s  Function 

In  each  of  the  procedures  to  be  discussed  in  this  chapter  it  is  necessary  to 
evaluate  an  integral  of  the  form 

l<(r')tST-f  )<n"  (8.1) 

where  g(r*-7' )  is  the  Green’s  function  described  in  Chapter  2  or 

-N0(R)+jJ0(R) 


kCm*  )  = 


(6.2) 


and  R  =  |  77'  | .  This  integral  in  equation  (6.1)  would  be  straightforward 
except  for  the  fact  that  N0(R)  tends  to  oo  as  R  goes  to  zero. 

The  effect  of  this  singularity  is  further  complicated  since  for  computer 
simulation  it  is  necessary  to  sample  the  function  at  discrete  positions.  If 
straightforward  sampling  is  performed  on  such  a  function  any  small  change  in 
the  location  of  the  sampling  grid  would  cause  a  large  change  in  the  sampled 
value  at  the  origin. 


Fortunately  it  is  possible  to  derive  an  approximate  value  for  the  sample  at 
the  origin.  The  actual  value  at  the  origin  is  not  so  important  but  when  the 
Green’s  function  is  multiplied  by  the  sampled  kernel  and  integrated  the  result 
should  be  identical  to  a  sampled  version  of  the  continuous  integral. 

Since  all  functions  have  been  sampled  they  can  be  assumed  to  be 
reasonably  smooth  and  the  integral  approximated 

Ef(ri)e(rj-7i)-  (M 

n 

When  7j  is  equal  to  7j  the  area  around  the  singularity  is  evaluated  and  the 
contribution  of  the  7j  term  of  the  summation  should  be  equal  to  the  original 
integral  in  a  square  region  around  the  origin. 

All  functions  are  sampled  with  interval  T,  in  both  the  x  and  the  y 
directions,  so  the  following  equality  is  necessary 

X  X 
2  2 

f^lgtfjTi)  =  /  /  r(f')g(77')df'.  (6.4) 

X  X 
2  2 

If  f(7j)  is  smooth  enough  then  it  can  be  considered  a  constant  within  the  small 
region  of  the  integral  and  brought  out  of  the  integral.  The  sample  of  the 
Green’s  function  is  now  written 


(6.5) 


2  2 


/  /  gfr?  W  ■ 
x_x 
2  2 


While  it  is  possible  to  evaluate  this  integral  using  approximations  to  the 
Bessel  functions,  a  better  solution  is  to  evaluate  the  integral  numerically.  The 
approximate  expansions  for  the  Bessel  functions  are  only  valid  for  small  values 
of  R  and  are  no  longer  accurate  as  the  sampling  interval,  T,  approaches  a 
quarter  of  a  wavelength. 


A  better  solution  is  to  numerically  integrate  this  function  by  dividing  the 
original  TxT  region  into  a  NxN  grid  of  points  and  summing  up  each  of  the 
samples.  The  samples  of  the  real  part  of  gfF-7' )  are  shown  in  Figure  6.1  for  a 
grid  size  of  4x4  and  16x16.  As  shown  in  this  figure  the  point  exactly  at  the 
origin  can  be  ignored  by  offsetting  the  grid  so  that  the  singularity  at  the  origin 
is  never  sampled.  The  effect  of  ignoring  this  undefined  location  can  then  be 
seen  by  examining  the  progression  of  average  values  as  the  number  of  samples 
in  the  T2  region  near  the  origin  is  increased.  These  values  are  shown  in  Table 


6.1  for  a  sampling  interval  of  —  wavelengths. 

4 


Notice  that  this  is  a  Cauchy 


sequence  and  converges  to  an  average  value  for  the  Creen’s  function  over  the 
TxT  region.  It  is  this  value  that  will  be  used  as  the  value  of  the  Green’s 
function  at  the  origin  in  each  of  the  computational  procedures  to  be  discussed 
in  the  remainder  of  this  chapter. 


6.3  Fixed  Point  Methods 

The  fixed  point  method  is  the  most  straightforward  computational 
approach  to  solve  an  equation.  In  the  remainder  of  this  section  it  will  be  used 
to  find  a  solution  of  both  the  Born  integral  and  the  Ricatti  equation. 

Consider  a  fixed  point  solution  to  the  equation  x  =  f(x)  defined  over  the 
region  x£[a,b],  Given  an  initial  guess,  x0  a  better  estimate  of  the  solution, 
can  be  found  by  [Sto80] 


Table  8.1.  Average  value  of  the  Green’s 


function  over  a  —X  region  sampled  on  an  NxN  grid. 

4 


N 

G(0) 

4 

0.0925259  +  jO.226659 

8 

0.0927197  +  jO.  225568 

16 

0.0927666  +  j0.225296 

32 

0.0927782  +  j0.225228 

64 

0.0927811  +  jO. 225211 

128 

0.0927818  +  jO.  225207 

256 

0.092782  +  jO. 225206 

512 

0.092782  +  j0.225206 

1024 

0.092782  +  jO. 225206 

x,  =  f(x0) 

(6.6) 

or  more  generally 

Xi  +  ,  =  f(Xj). 

(6.7) 

This  is  illustrated  in  Figure  6.2. 

While  this  method  doesn’t  always  converge,  it  is  possible  to  show  that  it 
will  converge  to  a  unique  solution  if  the  absolute  value  of  the  function’s  first 
derivative  is  always  less  than  one.  Mathematically  this  condition  can  be 
written 

|  r  (x)|  <1  xG[a,b],  (6.8) 

If  this  condition  is  true  then  it  is  also  possible  to  show  that  f()  is  a  contraction 
operator  or 

|  f(x,)-f(x2)|  <K|  X|~x2|  x£[a,b]  (6.9) 

where  K<  1. 

That  these  two  conditions  are  equivalent  is  easy  to  show  by  considering 
two  cases.  First,  assume  that  f()  is  a  contracting  operator.  From  the  definition 
of  the  derivative 


<  r 


l.iM  MJjJj  <,«  «, 


r  ♦,»  »j  i 


♦-<  t.«  U  Lt  »■»  «.•  «, 
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|f(x)|  =  lim 
'  1  h-o 


f(x4-h)-f(x) 


(x  +h)-x 

But  from  equation  (6.9)  this  ratio  and  the  derivative  are  less  than  K. 
To  show  the  converse  or  that 


(6.10) 


|  r(x)J  <1  xG[a,b)  (6.11) 

implies  a  contracting  operator  first  assume  the  opposite  case 
|  f(x2)-f(X|)|  >  |  x 2— x 1 1 .  Without  loss  of  generality  assume  x2  is  greater  than 
Xj  and  write  the  following  integral 

*2 

|  f(x2)-f(x ,)|  =  jF(x)dx.  (6.12) 

There  are  two  cases  to  consider.  First  assume  that  f(x2)~f(x,)>0  and  write 

*2 

0  <  f(x2)-f(x,)  =  /F(x)dx.  (6.13) 

*i 

But  if  F(x)<K  then  this  integral  can  be  bounded  from  above  by 

*2 

o  <  I  f(x2)-f(x,)|  =  / r  (x)dx  <  (x2-x,)K.  (6.14) 

*1 

This  contradicts  the  original  assumption. 

For  the  opposite  case,  f(x2)-f(xj)  <  0, 

*2 

|  f(x2)-f(x,)|  =  f(x,)-f(x2)  =  -fr  (x)dx.  (6.15) 

*i 

But  if  F(x)>-K  then 

*2 

0  >  f(x2)-f(x,)  =  JP(x)dx  >  -(x2-X|)K.  (6.16) 

*i 

For  both  cases,  f(x,)>f(x2)  and  f(x2)>f(xj),  the  original  assumption  is 
contradicted  therefore  proving  that  if  |F(x)|<K  then  the  function  is  a 
contracting  operator. 

To  study  the  convergence  of  the  fixed  point  method  to  the  correct 
solution,  £,  assume  a  value  for  Xj.  Then  the  error  after  the  i  +  1  iteration  is 
given  by 


Error  -  |  xi  +  r-Sj  =  |  f(x5)— . 


(6.17) 


But  since  £  is  a  fixed  point  solution  of  the  function  f  this  expression  for  the 
error  can  be  written 

Error  =  |  xi  +  1-£|  =  |  f(x;)-f(^)|  .  (6.18) 

By  the  contraction  property  of  Equation  (6.9)  this  last  term  is  bounded  from 
above  by 

Ixi+r£|  =  |f(x,)-f(0|  <K|xrf|.  (6.19) 

Thus  the  error  after  the  i  -E  1th  iteration  will  always  be  less  than  the  error  after 
i  iterations.  In  addition  since  K  is  less  than  one  the  sequence  j  x  -£|  is 
bounded  from  above  by  the  converging  geometric  series  x0K‘  therefore  the  fixed 
point  method  converges  to  a  correct  answer. 

That  this  solution  is  unique  can  be  seen  by  assuming  two  solutions,  and 
f2.  to  the  equation  x=f(x).  But  these  two  solutions  violate  the  contraction 
property  since 

I  «,-e2|  =  I  >  k)  e,-f2|  (6.20) 

Thus  if  f()  is  a  contracting  operator  then  the  fixed  point  solution  will  converge 
to  the  correct  answer. 

6.3.1  The  Born  Series 

Recall  from  Chapter  2  the  scattered  field  can  be  written 

ustH  = /gfr-T' )o(T' )u0(T' )df'  +  /gt?-r' Joff '  )us(r' )df ' .  (6.21) 

The  first  Born  approximation,  as  already  described,  simply  represents  the  first 
iteration  of  the  fixed  point  method  with  the  initial  guess  uj°)=0.  Thus  the  first 
iteration  is  written 

«s01  “  uB(T)  =  / gf?-?'  )o(T'  )u0(f'  )df ' .  (6.22) 

If  the  kernel,  gf?-?'  )o(r ' ),  is  a  contracting  operator  an  even  better  estimate 
can  be  found  by  substituting  u0(7)  +  uB(r)  for  u0(F)  in  equation  (6.22)  to  find 

^2|01  =  fg(?-7'  )o(T'  )[u0(F' )  +  uB(F'  ))df ' .  (6.23) 

In  general  the  i-th  order  Born  field  can  be  written 

=  /gO’-T'  WF'  )[u0(r' )  +  uffff  )|<tT' .  (6.24) 
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An  alternate  representation  is  possible  if  the  total  field  is  written 

ufH  “  u0(T)  +  Ui(r)  +  u2fF)+  •  •  •  (6.25) 


where 


Ufi  +  ijf7)  =  )o(T' )gfF-T' )df' . 


(6.26) 


By  expanding  equation  (6.24)  it  is  possible  to  see  that  the  fixed  point 
approximation  for  the  scattered  field,  u^,  is 


and  in  the  limit 


4’(n  =  E"iW 

j=0 


u(?)  =  u0(r)  +  u1(7)  +  u2(r)  +  u3(r)  +  •  •  •  . 


(6.27) 


(6.28) 


This  representation  (6.28)  has  a  more  intuitive  interpretation  analogous  to 
the  Huygens  principle.  The  Green’s  function  gives  the  scattered  field  due  to  a 
point  scatterer  and  thus  the  integral  of  equation  (6.24)  can  be  interpreted  as 
calculating  the  first  order  scattered  field  due  to  the  field  Uj.  For  this  reason  the 
first  order  Born  approximation  represents  the  first  order  scattered  field  and  U; 
represents  the  i’th  order  scattered  field. 

A  variation  of  the  integral  Born  series  presented  here  was  first  described 
and  implemented  by  Azimi  and  Kak  (Azi83).  In  [Azi83]  the  scattered  fields  are 
calculated  for  an  object  that  consists  of  multiple  cylinders  by  considering  the 
interactions  of  the  scattered  field  with  the  other  objects.  Using  the  equations 
presented  in  Chapter  5  and  (Wee64)  and  [Mor68]  it  is  possible  to  calculate  the 
exact  scattered  field  from  a  cylindrical  object  illuminated  with  a  plane  wave. 
Unfortunately  it  is  not  possible  to  calculate  the  exact  fields  scattered  by  a  pair 
of  cylinders  using  this  approach  because  the  field  from  one  object  interferes 
with  the  other. 

Instead  Azimi  and  Kak  propose  a  multiple  scattering  approach  where  the 
incident  field  is  first  scattered  against  each  cylinder  and  then  the  scattered 
fields  from  each  cylinder  are  propagated  to  the  other  cylinder(s)  where  they  are 
scattered  again.  This  is  illustrated  in  Figure  6.3  where  the  incident  field  is 
denoted  by  u0  and  the  field  Ujjk  denotes  the  field  that  has  scattered  off  of  object 
i,  then  object  j  and  finally  object  k. 

Since  the  field  scattered  by  a  cylinder  is  not  a  plane  wave  the  key  to  this 
procedure  is  to  calculate  the  scattered  fields  along  lines  between  the  cylinders. 
The  field  along  this  line  can  then  be  decomposed  into  plane  waves  [G0068]  and 
each  plane  wave  scattered  separately  by  the  cylinder. 
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While  nothing  is  known  about  the  convergence  of  this  series  its  use  in 
diffraction  tomography  is  limited  because  it  is  only  practical  when  the  field 
scattered  from  each  object  can  be  computed  exactly.  Thus  it  has  been  used  as 
a  method  to  generate  data  for  testing  diffraction  tomography  algorithms. 

Under  both  the  Dorn  and  Rytov  approximations  an  integral  of  the  form 

«i  +  if?)  =  JuifT'  )ofr '  )g(7-7'  )tfr'  (6.25) 

is  to  be  evaluated.  A  naive  approach  to  this  integral  would  be  to  evaluate  it 
numerically  but  doing  this  requires  on  the  order  of  N4  operations.  Performing 
the  integration  over  a  128  x  128  grid,  for  example,  would  require  over  270 
million  operations  per  iteration.  This  is  clearly  unreasonable. 

The  computational  requirements  can  be  greatly  reduced  by  noting  that  the 
Green’s  function,  g{7-7' ),  is  only  a  function  of  the  difference  between  the  two 
points  and  that  the  integral  can  also  be  interpreted  as  a  convolution.  Thus 
representing  the  convolution  in  the  frequency  domain  allows  for  an  efficient 
implementation  requiring  only  6N2logN  operations  to  do  the  integration.  For  a 
128  x  128  grid  this  is  only  700  thousand  operations  or  a  reduction  in  the 
computational  complexity  by  almost  400.  This  approach  is  efficient  using  an 
array  processor  and  for  an  object  sampled  over  a  64  x  64  grid  this  integral  can 
be  computed  in  under  2  seconds  using  a  Floating  Point  Systems  AP120B  array 
processor. 

To  evaluate  the  Born  integral  on  a  digital  computer  it  is  necessary  to  use  a 
truncated  and  sampled  version  of  the  object  function  and  the  field.  Truncating 
the  two  functions  is  error  free  because  the  object  function  is  assumed  to  have 
finite  support.  In  addition  since  only  a  finite  number  of  receivers  are  present 
this  further  limits  the  number  of  data  points  that  need  to  be  calculated.  On 
the  other  hand,  since  any  function  with  finite  support  has  infinite  bandwidth  it 
is  not  possible  to  sample  the  data  without  introducing  errors.  For  smoothly 
varying  data  it  is  possible  to  approximate  the  data  with  discrete  samples.  In 
the  work  to  follow  the  samples  are  taken  over  a  rectangular  grid  using  different 
sampling  intervals,  all  less  than  a  wavelength. 

If  the  series  Uj  decays  to  zero  then  the  total  field  is  given  by  a  summation 
of  each  scattered  field,  or 


Two  different  studies  were  performed  to  verify  this  approach  to  solving  the 
wave  equation.  Most  importantly  it  was  necessary  to  verify  that  the  total 
scattered  field  converged  to  the  same  answer  as  predicted  by  an  exact  solution 
to  the  wave  equation.  In  addition  it  was  necessary  to  determine  the  region  of 
convergence  of  the  Born  series.  These  issues  will  be  discussed  later  in  this 
chapter. 

As  mentioned  before,  the  integral  equation  in  (6.25)  is  efficiently  evaluated 
by  implementing  the  convolution  in  frequency  domain.  Recapitulating  this 
discussion,  the  frequency  domain  implementation  can  be  summarized  as 
follows. 


First  for  all  Xj  and  yk,  the  scattering  potential,  S,  is  calculated  from  the 
product  of  the  “incident”  field  and  the  object, 


s(Xj,yk)  =  ui(xj.yk)°(xj.yk)-  (6.30) 

Then  by  using  two-dimensional  FFT’s,  determine  the  following  two  transforms 


S(Uj,wk)  =  T2  FFT  (S(Xj,yk)) 


(6.31) 


Gi(uj,wk)=T2  FFT  (G(Xj,yk))  . 


(6.32) 


To  find  the  Fourier  transform  of  the  scattered  field  form  the  following  product 
in  the  frequency  domain  for  all  Uj  and  wk 

Ui  +  1(uj,wk)  =  S(uj,wk)  G(uj;wk).  (6.33) 


The  i  +  l’th  scattered  field  is  then  found  by  inverting  the  above  expression  to 
find 


i(xj.yk) 


IFFT 


'  ui+i(uj>wk)  ■ 


(6.34) 


To  properly  calculate  the  integral  using  the  FFT  approach  it  is  necessary 
to  remember  that  the  discrete  multiplication  implements  a  circular  convolution. 
A  circular  convolution  can  be  turned  into  an  aperiodic  convolution  by  zero 
padding  the  data  (Opp75j.  For  example  a  Floating  Point  Systems  (FPS) 
AP120B  Array  Processor  with  65,536  words  of  main  data  memory  can  only  deal 
with  arrays  up  to  128  x  128  elements.  This  memory  limitation  and  the  need  to 
implement  an  aperiodic  convolution  limits  the  evaluation  of  the  field  to  a  64  x 
64  grid. 
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The  exact  scattered  field  from  a  dielectric  cylinder  with  plane  wave 
illumination  is  well  known  in  the  literature.  This  exact  solution  to  the  wave 
equation  can  be  used  to  check  the  results  of  the  Born  iteration.  Figure  6.4 
shows  a  simulated  experiment.  In  this  experiment  a  plane  wave  is  incident  on 
a  cylinder  of  radius  2A  and  a  refractive  index  of  1.1.  The  scattered  field  was 
calculated  at  64  receiver  positions  along  a  line  7.75A  from  the  center  of  the 
object. 

Figure  6.5  show  the  exact  scattered  field  along  the  receiver  line.  In  each 
of  the  plots  that  follow  the  real  components  of  the  field  will  be  shown  with  a 
solid  line,  while  the  dotted  lines  represent  the  imaginary  part  of  the  field. 

Figure  6.6  shows  the  result  of  iterating  the  Born  integral.  The  first  Born 
approximation  (Figure  6.6a)  gives  a  very  poor  estimate  of  the  exact  field  since 
this  estimate  is  based  only  on  first  order  scattering. 

In  order  to  accurately  calculate  the  scattered  fields  it  is  necessary  to 
include  the  higher  order  scattered  fields.  This  is  shown  in  figure  6.6.  Clearly  in 
this  case  even  the  30’th  to  100’th  order  scattered  fields  are  significant  to  the 
total  field. 

By  comparing  figures  6.6  and  6.5  it  can  be  seen  that  in  this  case  the  Born 
series  converges  to  the  exact  scattered  fields.  This  simple  example  shows  the 
correctness  of  the  Born  iteration  code. 

The  Born  integral  defines  an  infinite  series  of  partial  scattered  fields  that 
are  summed  to  find  the  total  scattered  field.  An  important  measure  of  any 
infinite  series  is  its  region  of  convergence.  The  region  of  convergence  defines  a 
class  of  objects  where  the  Born  iteration  converges  to  the  exact  scattered  field. 
For  an  arbitrarily  complex  object  the  region  of  convergence  is  defined  over  an 
infinite  dimensional  space  since  an  infinite  number  of  parameters  are  needed  to 
describe  the  object. 

There  are  two  trivial  objects  that  can  be  analyzed  analytically.  First 
consider  a  pair  of  point  scatterers  separated  by  a  distance  of  R  located  at  7, 
and  72.  If  each  scatter  has  an  area  (or  volume)  of  a  and  the  object  function  for 
each  of  these  scatterers  is  equal  to  O  then  an  approximation  to  the  first  order 
Born  field  at72  due  to  the  scatterer  at  7,  is 

Uj*(72)  =s<tG(R)0  (6.35) 

Clearly  there  exists  a  value  of  O  such  that  U;1  will  have  a  magnitude  greater 
than  the  incident  field.  This  field,  U;\  can  then  be  scattered  with  the  point 
scatterer  at72  and  measure  the  field  at7j  to  find 
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Figure  6.4 


An  experiment  used  to  illustrate  the  higher  order  Born  series. 
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Figure  6.5 


Exact  scattered  field  for  the  object  shown  in  Figure  6.4.  The 
real  component  of  the  field  is  shown  as  a  solid  line  while  the 
imaginary  component  is  dashed. 
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Figure  6.6 


The  scattered  field  as  calculated  from  the  Born  series  using 
terms  numbered  from  1  to  100. 
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If  the  term  <rG(R)0  is  greater  than  one  then  each  succeeding  term  of  the  series 
Uj  will  diverge.  This  analysis  has  ignored  the  effect  of  the  other  scatterer  and 
the  field  at  T,  caused  by  the  scatterer  at  ?j  but  this  doesn’t  change  the  basic 
conclusion,  that  the  Born  series  can  converge  or  diverge  depending  on  the 
object. 


This  effect  can  also  be  analyzed  in  the  frequency  domain.  By  taking  the 
Fourier  transform  of  equation  (6.25)  the  convolution  with  the  Green’s  function 
can  be  expressed  as  a  multiplication  in  the  spatial  frequency  domain.  Equation 
(6.25)  becomes  [Sla83,  Sla84j 


ui+](R)  =  G(R) 


6|R).ui(K) 


(6.37) 


where  V  represents  convolution  in  the  frequency  domain. 

If  the  object  function  is  assumed  to  be  constant  for  all  space  then  O(R) 
becomes  an  impulse  and  equation  (6.37)  above  becomes 

Um(K)  =  OG(R)U,(K).  (6.38| 


This  is  a  simple  result  and  it  is  easy  to  see  that  the  series,  u^  will  diverge 
if  there  is  any  frequency  where 

G(R)0  >  1.  (6.30) 


In  this  simulation  study  the  region  of  convergence  for  a  single 
homogeneous  cylindrical  object  was  examined.  Since  any  cylinder  can  be 
completely  described  by  its  size  and  refractive  index  the  region  of  convergence 
is  defined  over  a  two  dimensional  space.  More  complicated  objects  could  be 
studied  but  the  results  would  not  be  as  graphical. 

It  has  already  been  shown  that  for  either  small  cylinders  or  small  changes 
in  the  refractive  index  the  first  Born  approximation  provides  a  good  estimate 
to  the  exact  scattered  field.  This  is  equivalent  to  saying  that  the  higher  order 
scattered  fields  decay  quickly  towards  zero.  If  is  easy  to  see  that  under  these 
two  conditions  the  Born  integral  series  will  quickly  converge  to  the  exact 
solution. 

The  total  energy  in  the  two  dimensional  field  is  used  as  a  simple  measure 
of  convergence.  Obviously  if  the  total  energy  in  the  i’th  scattered  field  is 
decaying  towards  zero  as  the  Born  integral  is  iterated  then  the  series  is 


converging.  On  the  other  hand  if  the  total  energy  is  increasing  then  the  Horn 
integral  can  not  possibly  converge. 

For  a  given  object  radius  the  region  of  convergence  is  determined  by 
conducting  a  binary  search  for  the  largest  refractive  index  where  the  Born 
series  converges.  For  each  combination  of  object  size  and  refractive  index  it  is 
necessary  to  make  a  decision  of  convergence  or  divergence  and  adjust  the 
refractive  index  accordingly. 


The  decision  of  convergence  or  divergence  is  made  by  studying  the  total 
energy  in  the  partial  scattered  field  as  the  iterations  are  performed.  For  this 
purpose  a  series  is  defined  to  be  convergent  if  during  each  of  four  iterations  the 
total  energy  in  the  partial  scattered  fields  is  monotonically  decreasing.  While  if 
the  energy  is  monotonically  increasing  then  it  is  decided  that  the  series 
diverges.  As  long  as  the  last  four  terms  are  not  monotonic  than  the  iterations 
continue.  The  energy  versus  iteration  number  for  the  experiment  of  figure  6.4 
is  shown  in  figure  6.7. 

Figure  6.8  displays  the  region  of  convergence  for  sampling  intervals  of 
0.125X,  0.25X  and  0.5X.  Each  plot  shows  the  maximum  refractive  index  as  a 
function  of  cylinder  radius.  For  all  experiments  with  a  refractive  index  below 
the  line  the  Born  series  converges,  while  for  all  experiments  above  the  line  the 
series  diverges. 

In  each  case  the  shape  of  the  curve  agrees  with  the  original  observation; 
the  Born  series  converges  for  either  small  objects  or  small  changes  in  the 
refractive  index.  The  dependence  of  the  region  of  convergence  on  the  sampling 
interval  is  still  under  study.  One  possible  explanation  is  that  the  numerical 
errors  are  larger  for  the  larger  sampling  intervals. 

Wrhile  more  complicated  objects  can  be  simulated  it  is  more  difficult  to 
present  the  region  of  convergence  in  a  simple  fashion.  In  general  the  region  of 
convergence  will  be  described  over  a  multidimensional  space  but  can  be 
reduced  to  two  dimensions  be  keeping  some  of  the  parameters  fixed.  Two 
simple  families  of  object  that  can  be  reduced  to  a  two  dimensional  space  will  be 
described  next.  In  each  case  the  results  will  be  compared  to  the  results  for  a 
cylindrical  object. 

A  simple  extension  of  the  previous  work  for  a  single  cylinder  is  to  consider 
two  cylinders  separated  by  a  small  distance.  Figure  6.9  shows  the  region  of 
convergence  for  two  orientations  of  the  cylinders  with  respect  to  the  incident 
field.  In  both  cases  two  identical  cylinders  are  separated  by  a  distance  of  1  X. 
Thus  the  region  of  convergence  is  reduced  to  two  dimensions,  the  radius  of  the 


Energy  in  Field 


two  cylinders  and  their  refractive  index.  The  solid  lines  in  Figure  6.9  show  the 
region  of  convergence  for  the  two  cylinders  compared  to  that  for  a  single 
cylinder  (shown  as  a  dashed  line). 

As  would  be  expected  the  region  of  convergence  for  two  cylinders  is  always 
smaller  than  that  for  a  single  cylinder.  This  is  true  since  the  convergence  is 
plotted  as  a  function  of  the  radius  of  one  of  the  cylinders  and  adding  another 
cylinder  increases  the  scattered  field.  As  already  described  increasing  the 
scattered  field  can  only  cause  the  Horn  series  to  converge  more  slowly. 

Figures  6.10  and  6.11  show  the  region  of  convergence  for  a  single  ellipse 
with  an  eccentricity  of  2  (length  of  major  axis  over  minor  axis  is  2.)  Since  an 
ellipse  can  not  be  described  by  its  radius  the  region  of  convergence  is  plotted  as 
a  function  of  the  length  of  both  the  major  and  minor  axis  of  the  ellipse.  The 
two  figures  differ  only  in  the  orientation  of  the  ellipse  with  the  incident  field 
and  in  both  cases  the  convergence  for  a  single  cylinder  as  a  function  of  its 
radius  is  plotted  as  a  dashed  line. 

The  convergence  of  the  Born  series  for  an  ellipse  (solid  lines)  compared  to 
a  cylinder  (dashed  lines)  is  consistent  with  the  idea  that  more  scattering  leads 
to  the  divergence  of  the  Born  series.  The  upper  solid  lines  show  the 
convergence  plotted  as  a  function  of  the  ellipse’s  major  axis  and  is  above  the 
dashed  lines  for  all  lengths.  This  is  because  an  ellipse  with  major  axis  of  length 
2r  has  less  area  than  a  circle  of  radius  r.  Conversely  the  lowest  line  (solid) 
plots  the  convergence  as  a  function  of  the  ellipse’s  minor  axis  and  is  always 
below  the  dashed  lines.  (Note  that  the  two  solid  lines  represent  the  same 
ellipse.  The  only  difference  is  that  region  of  convergence  is  plotted  with  respect 
to  different  parameters.) 

Figure  6.12  shows  the  same  data  as  the  previous  two  figures  but  now  the 
two  orientations  are  superimposed.  From  this  figure  it  is  easy  to  see  that  the 
convergence  of  the  Born  series  is  sensitive  to  the  orientation  of  the  ellipse. 
When  the  major  axis  of  the  ellipse  is  parallel  to  the  direction  of  the  incident 
field  the  region  of  convergence  is  reduced.  This  is  consistent  with  the 
limitation  described  in  Chapter  2  that  the  phase  change  in  the  field  as  it  travels 
through  the  object  is  a  good  indication  of  the  validity  of  the  Born 
approximation.  When  the  major  axis  is  aligned  with  the  incident  field  the 
phase  change  is  larger  than  when  it  is  perpendicular  to  the  field  and  thus  the 
Born  series  is  more  likely  to  diverge. 
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Figure  6.12  The  region  of  convergence  for  two  orientations  of  an  ellipse 

are  compared.  The  upper  curve  represents  the  convergence 
for  the  ellipse  shown  in  Figure  6.11  while  the  lower  curve 
represents  the  convergence  shown  in  Figure  6.10. 
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6.3.2  Born  Series  with  Attenuation 

While  to  this  point  only  fields  in  non-attenuating  media  have  been 
discussed  it  is  easy  to  also  talk  about  the  Born  series  when  attenuation  is 
present.  Now  the  refractive  index  and  the  wavenumber  are  no  longer  just  real 
valued  and  include  an  imaginary  component  to  represent  the  attenuation. 

The  relationship  between  the  real  and  imaginary  components  of  the 
wavenumber  are  easily  seen  by  examining  one  solution  to  the  homogeneous 
wave  equation 

u0(7)  =  eik°*  (6.40) 

where 7  =  (x,y).  If  k0  is  complex  and  equal  to 

k0  =  kr+jk;  (6.41) 

then  the  real  component  will  continue  to  represent  the  periodic  component  of 
the  field.  The  imaginary  component,  kj,  contributes  a  multiplicative  term, 

— kjc 

e  ,  that  causes  attenuation  of  the  plane  wave  with  increasing  distance. 

There  are  two  approaches  that  can  be  used  to  deal  with  attenuation.  In 
the  simpler  approach  the  average  wavenumber,  k0,  is  real  and  all  of  the 
attenuation  is  a  perturbation  from  the  average  wavenumber.  Thus  the  object 
function  is  complex  and  as  will  be  shown  shortly  the  region  of  convergence  is 
reduced  for  large  attenuations.  In  the  second  approach  the  attenuation  is 
included  in  the  average  wavenumber  thus  reducing  the  magnitude  of  the  object 
function.  The  only  difference  in  the  formulation  is  that  the  Green’s  function 
changes  but  now  increasing  attenuation  leads  to  a  larger  region  of  convergence. 

When  the  attenuation  of  the  object  is  treated  as  just  a  perturbation  of  the 
object  function  from  the  real  valued  wavenumber  then  the  effect  is  to  reduce 
the  region  of  convergence.  This  is  shown  in  Figure  6.13  where  it  can  be  seen 
that  the  magnitude  of  the  object  function,  not  just  the  refractive  index, 
determines  the  convergence  of  the  Born. 

The  work  first  done  with  the  Bom  approximation  assumed  that  the 
average  refractive  index  was  real  valued  only.  Since  the  real  part  of  the 
wavenumber  represents  the  speed  of  the  wave  and  the  imaginary  part  its 
attenuation,  any  attenuation  in  the  object  is  included  in  the  unknown 
perturbation.  The  magnitude  of  the  perturbation  determines  the  applicability 
of  the  Bom  and  Rytov  series  therefore  a  more  accurate  estimate  of  the  average 
refractive  index  will  lead  to  smaller  perturbations  and  better  results  with  first 
order  diffraction  tomography  algorithms. 
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The  region  of  convergence  for  the  Born  series  as  a  function  of 
cylinder  radius  and  attenuation.  The  attenuation  is  plotted 
in  units  of  nepers  (1  neper  represents  an  attenuation  of  the 
fields  amplitude  by  63%  per  wavelength.). 


Figure  6.13 


Consider  an  experiment  where  an  object  and  the  surrounding  media  are 
both  highly  attenuating.  This  might  be  typical  of  a  microwave  tomography 
experiment  where  the  attenuation  of  microwaves  is  predominately  due  to  the 
water  molecule. 

In  this  case  a  small  perturbation  model  would  be  more  accurate  if  the 
average  value  of  the  wavenumber,  k0,  is  assumed  to  include  an  imaginary 
component.  Thus  the  real  part  continues  to  represent  the  spatial  frequency  of 
the  wave  while  the  imaginary  part  indicates  the  bulk  attenuation  of  the  wave 
as  it  travels  through  the  media. 

An  important  part  of  this  derivation  is  to  remember  that  a  field  can  be 
described  in  two  different  manners;  if  you  like,  there  are  two  sets  of  basis 
functions  that  can  be  used.  A  field  that  satisfies  the  Helmholtz  equation  is 
usually  described  in  terms  of  plane  waves.  A  plane  wave  is  an  exponential 
solution  to  the  Helmholtz  equation  and  for  a  plane  field  described  by 

u(F)  =  e*'7  (6.42) 


where 


k  =  (kx,ky)  (6.43) 

then  valid  plane  waves  satisfy 

k2  +  k2  =  k02  (6.44) 

k0  in  this  equation  represents  the  wavenumber  of  the  media. 

A  problem  with  this  approach  is  that  both  kx  and  ky  can  be  complex.  The 
complex  valued  wavevectors  lead  to  evanescent  waves  which  attenuate  with 
distance.  While  for  most  applications  the  evanescent  fields  can  be  ignored 
(they  tend  to  be  small  compared  to  the  non  attenuating  components)  the  same 
assumption  can  not  be  made  when  calculating  the  field  inside  the  object  or 
when  the  object  is  in  an  attenuating  media. 

Consider  an  attenuating  plane  wave  propagating  in  the  x  direction 

upf)  =  etfr  =  ^  =  ei(Q+^x  /?>0  (6.45) 

In  this  case  a  represents  the  phase  term,  while  /?  represents  the  attenuation. 
For  positive  x  this  is  an  attenuating  plane  wave,  but  for  negative  x  the  wave 
grows  exponentially.  Thus  it  is  normally  necessary  to  specify  that  the  wave  is 
zero  for  x  <0. 

Due  to  the  non  symmetry  of  attenuating  plane  waves  and  the  efficiency  of 
Fast  Fourier  Transform  algorithm  a  much  more  natural  set  of  basis  functions  is 
provided  by  the  Fourier  domain.  In  this  approach  the  field  is  represented  as  a 


sum  of  Fourier  components.  If  the  field  does  not  have  any  evanescent 
components  the  Fourier  and  the  plane  wave  representations  are  identical. 

The  distinction  becomes  important  when  attenuating  fields  are  considered. 
Mathematically  the  two  approaches  are  equally  valid  but  while  an  attenuating 
plane  wave  has  a  single  component  in  the  plane  wave  representation  it  has  an 
infinite  number  of  Fourier  components.  The  difference  is  further  illuminated  if 
an  attenuating  plane  wave  is  propagating  through  a  homogeneous  media.  This 
field  is  represented  as  a  single  plane  wave  that  satisfies  the  wave  equation  but 
its  Fourier  representation  has  an  infinite  number  of  components.  While  each  of 
its  Fourier  components  do  not  satisfy  the  wave  equation  they  do  represent  a  set 
of  basis  functions  for  describing  linear  operators.  For  this  reason  the  Fourier 
approach  is  optimum  for  propagation  problems  [G0068]  and  convolution 
integrals. 

Again,  like  was  done  for  the  non  attenuating  Born,  Figure  6.14  shows  the 
components  of  the  Born  integral  for  a  complex  wavenumber.  The  major 
differences  are  that  the  incident  field  becomes  an  attenuating  plane  wave  and 
the  attenuating  Green’s  function  is  more  straightforward  to  calculate  since  it 
no  longer  has  a  singularity  in  the  frequency  domain. 

The  incident  field 


u(?)  =  x>0  (6.46) 

is  a  complex  (2  dimensional)  function  and  its  Fourier  transform  can  be  found 
by  considering  it  as  a  multiplication  of  a  complex  sinusoid  by  a  one  sided 
exponential.  The  following  one  dimensional  Fourier  pairs  are  used: 

e^Wot  <->  2^w-w0)  (6.47) 

and 


e-at  t>0  ~  1  ~ — .  (6.48) 

jw  +  a 

Multiplication  in  the  space  domain  corresponds  to  convolution  in  the  frequency 
domain  so  the  Fourier  transform  of  an  attenuating  sinusoid  is  written 


j(ur-w0)  +  a  ' 


(6.49) 


In  two  dimensions  the  Fourier  transform  of  the  incident  field  is  written 
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where  the  wavevector  of  the  incident  field  is  given  by  A  =  (a,/?).  This  is  shown 
in  Figure  6.14(Incident  Field). 

The  Green’s  function  for  complex  values  of  k0  is  even  simpler  than  the 
non-attenuating  case.  By  taking  the  Fourier  transform  of  the  wave  equation 
with  a  delta  forcing  function  it  is  easy  to  see  that 


Since  |  A |  2  is  real  and  k©  is  complex  there  are  no  singularities  in  this  function. 
The  Fourier  transform  of  the  Green’s  function  is  still  circularly  symmetric  and 
is  shown  in  Figure  6.14(Green’s  Function). 

The  Fourier  transform  of  the  two  dimensional  scattered  field  can  be 
written  now  as 

U,(R)  =  G(R)  I  U0(R)  .  0{R) ).  (6.52) 

The  convolution  of  the  incident  field  and  the  object  is  a  shifted  and  a 
smeared  version  of  the  Fourier  transform  of  the  object.  The  convolution  can 
be  considered  in  two  parts.  The  smearing  caused  by  the  width  of  the  incident 
field  in  the  spatial  frequency  domain  can  be  ignored  since  it  just  redistributes 
some  of  the  energy  of  the  object  function’s  Fourier  transform.  The  remaining 
component,  the  shift  in  the  frequency  domain,  is  identical  to  the  non 
attenuating  Born  case. 

As  derived  above,  the  Green’s  function  for  an  attenuating  media  does  not 
have  a  ring  of  singularities  and  therefore  samples  a  semi-circular  region  of  the 
modified  object  function.  This  is  shown  in  Figure  6.14(Scattered  Field). 

This  procedure  naturally  leads  to  the  Bom  series  for  objects  with  bulk 
attenuation.  Figure  6.15  shows  a  composite  graph  of  the  region  of  convergence 
for  a  number  of  different  attenuating  media  between  0  and  1  nepers  per 
wavelength.  It  is  interesting  to  note  that  the  region  of  convergence  gets  larger 
as  the  attenuation  increases.  This  is  due  to  the  reduction  in  multiple  scattering 
because  to  the  attenuating  term  in  each  wave.  Thus  for  an  average 
attenuation  of  1  neper  per  wavelength  (the  amplitude  of  the  field  drops  by  3dB 
per  wavelength)  the  Bom  series  converges  for  all  objects  with  a  refractive 
index  less  than  20%.  In  addition  the  convergence  of  the  Bom  is  less  sensitive 
to  the  size  of  the  object  since  the  waves  are  attenuated  before  they  travel  the 
complete  distance  of  the  object. 


0.3.3  Rytov  Series 

The  implementation  of  the  Rytov  series  is  much  like  that  of  the  Born. 
Thus  from  Chapter  2  the  Rytov  integral  for  the  scattering  phase  is  written 

*,tf)  =  -~l  s(M")  |v^(r')  v^,fr')+o(r')]  «.(i ")  Jr'.  (6.53) 

u0(7)  v 

A  fixed  point  solution  to  this  equation  is  possible  if  a  guess  for  is  used  in 
the  right  side  of  this  equation  and  a  new  value  for  the  scattered  phase,  is 

computed.  If  the  kernel  of  equation  (6.53)  is  a  contracting  operator  then 
will  be  a  better  estimate  of  the  scattered  phase.  This  iteration  step  can  be 
carried  out  as  often  as  desired  until  the  change  in  the  scattered  phase  is  small. 

The  computer  implementation  of  the  Rytov  series  is  more  difficult  than 
the  Born  due  to  the  derivatives  inside  the  integral.  In  the  Born  approximation 
it  is  possible  to  decompose  the  series  and  think  of  the  entire  iteration  as 
modeling  higher  order  scattering.  On  the  other  hand  in  the  Rytov  integral  the 
[V^jfP'  )  V^,(T' )]  term  is  a  non-linear  function  of  the  scattered  phase  and  thus 
it  is  not  possible  to  consider  each  term  separately. 

The  computer  implementation  of  the  Rytov  series  is  made  especially 
difficult  since  the  scattered  potential  does  not  have  finite  support  as  it  does  in 
the  Born  series.  Recall,  in  the  Born  series  the  scattering  potential  is  a  function 
of  the  product  of  the  scattered  field  and  the  object  function  and  the 
convolution  integral  need  be  evaluated  only  where  the  object  function  is  non¬ 
zero.  Thus  carrying  out  the  Born  integration  over  a  finite  region  does  not 
introduce  any  errors. 

This  is  not  the  case  with  the  Rytov  integral  since  the  scattering  potential 
is  now  given  by  the  expression 

lV^(?)  V^j(f  )]u0f?).  (6.54) 

In  general,  the  scattered  phase,  rf) ,,  is  not  equal  to  zero  so  limiting  the 
integration  to  a  finite  region  will  always  introduce  errors. 

Since  the  derivative  is  a  linear  operator  it  can  be  implemented  as  a 
convolution  integral.  Unfortunately,  as  will  be  shown  shortly,  the  structure  of 
the  problem  is  such  that  the  most  accurate  method,  based  on  FFT’s,  is  not 
workable. 

When  using  an  FFT  to  implement  a  convolution  integral  it  is  necessary  to 
zero  pad  the  original  data  so  that  the  FFT  represents  an  aperiodic  convolution. 
While  this  technique  works  very  well  for  most  signals  it  has  disastrous 
consequences  when  calculating  the  derivative  of  the  field.  Since  the  field  never 
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goes  completely  to  zero  there  is  always  a  sharp  transition  between  the  held  and 
the  start  of  the  zero  padding.  This  transition  leads  to  a  large  value  of  the 
derivative  at  this  point  and  eventually  to  large  errors  in  the  scattered  field. 

The  standard  procedure  for  dealing  with  this  problem  is  to  use  a  window 
to  smooth  out  the  transition.  This  solution  is  not  viable  here  since  the  problem 
is  to  calculate  the  scattered  phase  at  the  outer  edge  of  the  grid,  exactly  where 
the  effect  of  the  window  is  greatest.  Thus  even  though  the  kernel  for  the 
derivative  operator  is  very  compact  the  long  tails  lead  to  errors  with  an  FFT 
based  implementation. 

A  better  solution  is  to  approximate  the  derivative  operator  with  a  two 
point  kernel  and  make  appropriate  adjustments  at  the  edges  of  the  grid.  A 
third  order  polynomial  is  fit  to  the  three  points  (— t,y_  j),  (0,y0)  and  (t,y j )  with 
the  function  [Sto80] 


,  y-ix(x-t)  y0(x-t)(x  +  t)  y,(x  +t)x 
f,x)  =  ^ 

The  first  derivative  of  this  polynomial  is  found  to  be 

P(x)  = 
v  ’  2t 


(6.55) 


(6.56) 


At  the  edge  of  the  grid  all  three  values  of  the  field  are  not  available  and  so  a 
second  order  polynomial  is  fit  to  the  two  points  and  the  derivative  becomes 


yo-y-i 


(6.57) 


yryo 


(6.58) 


These  operators  are  straightforward  and  allow  the  V2  operator  to  be  computed 
quickly  in  the  time  domain. 

The  convergence  of  the  Rytov  series  is  much  like  that  of  the  Born.  That 
the  Rytov  series  converges  to  the  correct  answer  is  shown  in  Figures  6.16  and 
6.17.  Figure  6.16  shows  the  exact  scattered  field  from  a  2X  cylinder  with  a 
refractive  index  of  1.13.  The  field  was  measured  at  a  receiver  line  7.75X  from 
the  center  of  the  cylinder  and  sampled  every  1/4X.  Finally  Figure  6.17  shows 
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that  the  Rytov  series  does  converge  to  the  exact  solution  for  the  scattered  field. 

The  region  of  convergence  of  the  Rytov  series  is  compared  to  that  of  the 
Born  in  Figure  6.18.  A  number  of  works  [Kel69,  San70,  Man70  and  Wes84] 
have  discussed  this  issue  but  it  is  not  clear  from  the  theoretical  work  which 
series  is  superior.  These  numerical  implementations  of  the  Born  and  Rytov 
series  show  the  Rytov  to  converge  for  a  large  range  of  objects  than  the  Born 
does.  This  is  especially  surprising  considering  the  different  domains  of  validity 
of  the  Born  and  Rytov  approximations. 

The  difference  between  the  Born  and  Rytov  series  is  highlighted  in  Figures 
6.19  and  6.20.  These  two  figures  show  the  convergence  of  the  Rytov  series  for 
an  object  made  up  of  two  cylinders  and  an  elliptical  object.  While  in  the  Born 
approximation  the  orientation  of  the  object  changed  the  convergence  of  the 
series  the  same  is  not  true  for  the  Rytov  series.  The  convergence  results  for 
the  Rytov  series  are  identical  for  either  orientation  (0  and  90  degrees). 

The  behavior  of  the  Rytov  series  with  an  attenuating  object  is  shown  in 
Figure  6.21.  Like  the  Born  series  the  Rytov  series  is  relatively  insensitive  to 
attenuation  in  the  object  until  the  attenuation  becomes  large  enough.  The 
attenuation  at  which  the  convergence  of  the  Rytov  series  falls  to  zero  is 
dependent  on  the  radius  of  the  cylinder. 

If  the  object  and  the  media  have  an  average  attenuation  then  the  Rytov 
series  will  converge  more  easily.  This  is  shown  in  Figure  6.22.  Using  an 
attenuating  Green’s  function  reduces  the  field  at  distances  far  from  the  object 
and  thus  makes  it  easier  for  the  Rytov  series  to  converge. 

6.4  Matrix  Formulation 

An  alternative  to  the  fixed  point  methods  like  the  Born  and  the  Rytov 
was  shown  by  Kaczmarz  [Kac37],  applied  to  the  forward  scattering  problem  by 
Richmond  [Ric65]  and  extended  to  inverse  scattering  by  Johnson  [Joh83,  Tra83 
and  Joh84].  The  Kaczmarz  algorithm  has  found  widespread  use  in 
tomographic  imaging  based  on  ray  tracing.  Its  use  for  this  type  of  problem  and 
a  discussion  of  several  possible  optimizations  and  tricks  is  discussed  in  [Her73, 
Her76  and  Her80j.  While  the  Born  and  the  Rytov  series  use  discrete  math  to 
implement  a  continuous  solution  to  the  Helmholtz  equation,  a  different 
approach  is  possible  if  the  field  and  the  object  are  first  discretized.  The 
Helmholtz  equation  now  becomes  a  matrix  equation  and  with  appropriate 
manipulations  can  be  put  in  the  form 
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Figure  6.20  The  region  of  convergence  of  the  Rytov  series  for  an  ellipse. 

The  upper  solid  line  is  plotted  as  a  function  of  the  major  axis 
while  the  lower  solid  line  is  a  function  of  the  minor  axis. 
The  solid  line  represents  the  convergence  for  a  single 
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Figure  6.21  The  convergence  of  the  Rytov  series  is  shown  as  a  function 

of  the  cylinders  radius  and  attenuation.  The  attenuation  of 
the  object  is  shown  in  nepers. 
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Figure  6.22 


The  convergence  of  the  Rytov  series  is  shown  as  a  function 
of  the  cylinder  radius  and  the  average  attenuation  of  the 
media.  The  attenuation  is  shown  in  nepers. 
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where  x  is  the  continuous  field,  A  is  a  function  of  the  object  and  the  Green’s 
function  and  b  represents  the  effects  of  the  incident  field. 

The  Kaczmarz  algorithm  is  used  to  solve  this  matrix  equation  because  it 
operates  on  only  one  row  of  the  matrix  at  a  time.  This  property  is  important 
due  to  the  large  size  of  the  matrix  involved  (as  many  as  4000  equations  and 
unknowns.)  The  Kaczmarz  algorithm  belongs  to  a  class  of  operators  known  as 
Row  Action  Methods  and  is  described  by  Tanabe  [Tan7l]  and  by  Censor 
(Cen81). 

Most  work  with  large  matrices  assumes  that  a  significant  fraction  of  the 
matrix  is  zero.  Thus  it  is  only  necessary  to  store  the  location  and  the  value  of 
the  non  zero  elements  and  then  it  is  possible  to  use  sparse  matrix  techniques  to 
solve  the  equation  (see,  for  example,  Chapter  3  of  [Hey85].)  Since  every  element 
of  the  Green’s  matrix  is  non  zero  it  is  not  possible  to  use  these  techniques  and 
instead  the  Kaczmarz  approach  allows  the  calculation  of  a  solution  using 
modest  amounts  of  computer  memory. 

As  will  be  shown  later  the  Kaczmarz  algorithm  always  converges  to  a 
proper  solution  of  the  discrete  equation  Ax  =  b.  If  the  discrete  representation 
of  the  field,  the  object  and  the  Green’s  function  accurately  model  the  true 
functions  then  the  Kaczmarz  solution  will  satisfy  the  Helmholtz  equation. 

An  exact  solution  to  the  wave  equation  is  given  by  the  integral  equation 

us(7)  =  /off' )  (u0(?' )  +  us(f ' )]  g(?-7'  )df 1 ' .  (6.60) 

By  sampling  faster  than  the  Nyquist  rate,  each  of  the  terms  in  the  above 
equation  can  be  discretized  without  errors  as 

us1,J  =  us(iT,jT)  (6.61) 


u,j’i  =  uc(iT,jT) 
o*’j  =  o(iT,jT) 
gi’i  =  g(iT,jT). 

The  discrete  version  of  the  Helmholtz  equation  can  now  be  written 

«,u  =  E  E  «kJ  I  “o'1  +  «,kJ  ]  f  k'H 

k  =  -oo  I  =  -  oo 


(6.62) 

(6.63) 

(6.64) 


(6.65) 


Again,  as  in  the  implementation  of  the  Bom  and  Rytov  series,  assume  that  the 
object  has  finite  support.  Since  the  object  multiplies  each  term  in  the 


summation  the  summation  need  be  carried  out  only  for  those  values  of  k  and  1 
where  ok’’  is  non-zero.  Thus  without  loss  of  generality  the  object  will  be 
assumed  to  exist  only  for  l<i<N  and  l<j<N  and  all  summations  will  be 
assumed  to  go  from  1  to  N. 

The  forward  process  will  be  described  first.  To  do  this  rearrange  the 
discrete  version  of  the  wave  equation  to  find 

£E  ok-'  u0k-‘  g i"k-j-1  =  uj-i  -  ry;  ok'‘  uk-‘  gi  k'j'*.  (6.66) 

k  I  k  I 

The  Ug’j  term  can  then  be  moved  into  the  summation  to  find 

EE  °k’'  “o’1  6i“k,j"1  =  EE  (  #  -  1  uk-‘  (6.67) 

k  1  k  1 

where  the  Kroenker  delta,  is  defined  as 

1  for  i  =  k  and  j=l 

0  elsewhere  (6.68) 


To  put  this  in  standard  matrix  form  represent  each  field  as  a  one¬ 
dimensional  vector  as  follows 


x«  =  lx00’0,  x^’1, 


•  x°N -1  xi°  xP  •  •  ■  x‘-N  l 
>A0  >  A0  y  A0  »  >A0  y 


(6.69) 
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(6.70) 


yn-i,o  .  .  .  „N-|,N-1 

A  >  >A 

where  (]T  represents  vector  transpose.  An  N^xN2  matrix,  A,  can  now  be 
defined  that  represents  the  effects  of  all  summations.  Since  the  object  and  the 
Green’s  function  are  both  known  they  can  be  combined  into  a  single  matrix 
and  the  discrete  version  of  the  wave  equation  can  be  written 


Ax  =  b 

where  the  terms  of  the  A  matrix  are  given  by 


(6.71) 


A  = 


(6.72) 
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The  field,  x,  is  a  one  dimensional  matrix  with  N2  elements  and  the  constant 
vector  b  is  given  by 

b  =  Ax0  (6.73) 

At  first  glance  this  equation  (6.71)  represents  an  especially  simple  form  for 
the  scattered  field;  that  is  until  the  size  of  the  vectors  are  considered.  For  a 
small  64x64  reconstruction  the  b  and  x  vectors  have  4096  complex  elements 
and  the  A  matrix  is  a  square  4096x4096  matrix  with  more  than  16  million 
complex  elements.  Inverting  a  matrix  of  this  size  would  require  over  32  Mega- 
words  of  memory  or  more  than  that  which  exists  on  all  but  a  handful  of 
processors  today. 

There  are  two  tricks  to  solving  this  problem.  The  first  of  which  is  to 
realize  that  it  is  not  necessary  to  find  the  inverse  of  the  A  matrix  but  only  to 
find  a  vector  field,  x,  which  satisfies  the  discrete  wave  equation.  Secondly  by 
using  a  row  action  method  such  as  that  proposed  by  Kaczmarz  and  Johnson  et 
al,  it  is  no  longer  necessary  to  store  the  entire  A  matrix  in  memory.  Thus  it  is 
possible  to  solve  the  system  of  equations  storing  only  4N2  equations  at  a  time. 
For  a  64x64  field  this  represents  only  16000  elements  so  the  storage 
requirements  are  reduced  by  a  factor  of  1000. 

As  described  in  [Ros82],  the  Kaczmarz  method  finds  a  vector  x  that 
satisfies  the  equation  Ax  =b  by  considering  each  row  of  the  matrix  A  to 
represent  a  separate  equation.  Thus  the  matrix 

all  a12  a13 

A  —  a^  &22  ^3  (6.74) 

a31  a32  a33 

can  be  considered  to  be  three  equations  of  the  form 


(6.75) 


bl  =  allxl  +  a12x2  +  a13x3 
b2  =  ^21*1  a22x2  a23x3 
b3  =  a31x  1  +  a32x2  a33x3 


In  terms  of  an  n-dimensional  space  each  of  these  equations  represents  a  single 
hyperplane  and  the  intersection  of  all  ‘n’  planes  describes  a  single  solution 
point  x. 

The  Kaczmarz  algorithm  iteratively  refines  an  initial  guess  by  projecting 
the  point  onto  each  hyper-plane  in  sequence.  The  process  of  computing  the 
projection  onto  the  hyperplane  also  represents  finding  the  point  on  the 
hyperplane  closest  to  the  original  guess.  As  will  be  shown  later  this  new  x  will 
always  lie  closer  to  the  new  solution  vector  than  the  original  guess. 

If  the  i’th  row  of  the  A-matrix  is  denoted  as  a;  and  <  ,  >  represents  the 
dot  product  then  a  better  solution  to  the  equation  AxJ=:b  is  given  by 


x 


j  +  i 


<a5)xi>-bi 
xJ - a;. 

<ai,ai> 


(6.76) 


This  operation  is  illustrated  in  Figure  6.23. 

Simple  geometrical  arguments  should  convince  the  reader  that  this 
equation,  (6.75)  will  always  produce  a  better  estimate  of  the  solution  vector  x. 
For  a  two  dimensional  case  this  situation  is  illustrated  in  Figure  6.24.  In  this 
example  a  point,  xlt  on  the  line  CD  will  first  be  projected  onto  the  line  AB  or  a 
solution  of  the  equation 


aHx  1  +  ai2x2  “  bl  (6.77) 

For  any  point  on  the  line  CD,  the  point  x2  on  the  line  AB  is  always  a  better 
estimate  of  the  solution  then  the  original  point.  Thus  the  Kaczmarz  algorithm 
always  converges.  The  distance  between  two  points  Xj  and  Xj  +  ,  will  be  defined 
in  the  normal  Euclidean  sense  or 

Distance2  =  <xj-x^  +  \x^-x'  +  1>.  (6.78) 

Since  the  solution  vector,  x,  lies  along  the  line  AB  an  initial  estimate  can 
always  be  improved  by  projecting  the  point  onto  the  line.  The  point  x^*1  is 
closer  to  the  solution  than  all  other  points  on  the  line  CD  so  this  always 
represents  a  better  estimate.  Since  all  points  on  the  line  CD  project  to  the 
point  xJ  +  1  this  procedure  always  reduces  the  error.  Thus  it  is  easy  to  see  that 
the  error  is  monotonioly  decreasing  and  in  addition  will  always  converge  to 
zero.  (For  now  the  cases  where  the  system  is  either  over  determined  or 
underdetermined  are  ignored.  Both  of  these  cases  represent  systems  of 


equations  that  don’t  necessarily  have  a  unique  solution.) 

The  speed  of  convergence  is  proportional  to  the  independence  of  the  rows 
of  the  A  matrix.  Figure  6.25  shows  the  convergence  for  two  widely  separated 
cases.  In  Figure  6.25a  the  two  hyperplanes  are  perpendicular  and  the 
Kaczmarz  algorithm  converges  to  the  solution  in  one  iteration,  no  matter  what 
the  starting  point  is.  (One  iteration  is  defined  as  projecting  the  x  vector  onto 
each  of  the  hyperplanes.)  On  the  other  hand  in  Figure  6.25b  the  hyperplanes 
are  nearly  parallel  and  it  will  take  Kaczmarz  algorithm  many  iterations  to 
converge  to  the  correct  solution. 

Since  the  Kaczmarz  method  only  works  with  a  single  row  of  the  A  matrix 
at  a  time  it  is  possible  to  make  a  space-time  tradeoff.  The  form  of  the  A 
matrix  is  simple  enough  that  it  is  relatively  inexpensive  to  recompute  each  row 
of  the  matrix  as  it  is  needed.  While  on  many  computer  systems  it  is  possible 
to  precompute  the  A  matrix  and  store  it  on  disk,  getting  the  data  back  off  can 
be  time  consuming.  For  the  computers  accessible  at  Purdue  (Floating  Point 
Systems  AP120B  and  Control  Data  Corporation  Cyber  205)  it  is  more  cost 
effective  to  recompute  the  A  matrix  as  needed.  This  approach  has  made  the 
problem  solvable  for  more  than  trivial  sized  matrices. 


For  an  implementation  of  the  forward  scattering  problem  on  a  CDC  Cyber 
205  super  computer  calculating  one  row  of  the  A  matrix  requires  98,000 
floating  point  operations  (64x64  grid).  Since  calculating  the  projection  requires 
57,000  floating  point  operations  this  implementation  takes  2.7  times  longer 
than  the  ideal  (all  values  of  the  A  matrix  available  immediately.)  On  the  other 
hand  retrieving  the  data  off  disk,  either  explicitly  or  using  virtual  memory, 
could  take  100  times  longer  than  the  ideal  situation*.  Thus  for  this  system  of 
equations  the  tradeoff  is  easy  to  make.  With  a  Cyber  205  one  iteration  of  the 
Kaczmarz  algorithm  takes  233  million  floating  point  operations  and  can  be 
computed  in  1  second  of  CPU  time.  This  represents  a  real  cost  of  21*  per 
iteration  as  charged  at  the  standard  rate  by  the  Purdue  University  Computer 
Center. 

To  make  the  implementation  of  the  Kaczmarz  algorithm  as  fast  as  possible 
several  quantities  are  precomputed  and  stored  for  quick  access.  The  x  and  y 
coordinates  of  each  grid  point  do  not  change  during  the  course  of  the  problem 


Theoretically  it  should  be  possible  to  organize  the  data  on  the  disk  in  a  sequential 
fashion  and  overlap  execution  time  with  IO  time.  With  high  enough  bandwidth  the  total 
execution  time  should  be  nearly  equal  to  the  ideal  situation.  Unfortunately,  most  operat¬ 
ing  systems  are  not  this  cooperative. 


The  orthogonality  of  the  hyperplanes  determines  the  rate  of 
convergence.  If  the  hyperplanes  are  perpendicular  then  the 
solution  will  be  reached  in  only  one  iteration  (a)  while  it  will 
take  much  longer  if  the  hyperplanes  are  nearly  parallel  (b). 


.*  y  ^  *-«  ja  is  >^  tL  i*»  **.' >« .  tC A,nL*ti.4ti,^i.1' 


■wwwuuumKvuvv«*w 


162 


and  thus  it  is  possible  to  store  these  in  main  memory  and  not  recalculate  them. 
In  addition  the  Green’s  function  is  circularly  symmetric  and  thus  a  function  of 
only  the  radial  distance.  By  precomputing  the  values  of  the  Green’s  function 
along  a  single  radial  it  is  possible  to  use  bilinear  interpolation  to  quickly 
compute  the  Green’s  function  for  each  grid  point  given  the  radial  distance  and 
the  values  along  the  radial. 

The  algorithm  for  one  iteration  of  the  Kaczmarz  algorithm  can  be  written 
as 

For  each  equation  (representing  the  scattered  field  at  a  single 
point) 

Compute  the  Green’s  Function 

•  Subtract  a  vector  representing  the  x  position  of  each 
grid  point  from  the  current  equation. 

•  Subtract  a  vector  representing  the  y  position  of  each 
grid  point  from  the  current  equation. 

•  Square  each  distance  and  add 

•  Find  square  root  and  multiply  by  scaling  factor  to  find 
the  argument  of  the  Green’s  function. 

•  Use  bilinear  interpolation  to  find  the  Green’s  function 
at  every  point. 

Find  the  A  matrix 

•  Multiply  the  Green’s  function  and  the  object 

•  Subtract  from  the  identity  matrix  (fy) 

Project  xJ  onto  the  hyperplane. 

An  additional  complication  in  this  approach  is  caused  by  the  use  of 
complex  numbers.  While  the  dot  product  operation  is  defined  for  complex 
vectors  better  results  are  obtained  if  each  complex  equation  is  considered  to  be 
two  real  valued  equations.  This  simple  change  reduces  the  error  by  a  factor  of 
100  or  more  but  does  increase  the  number  of  equations  for  a  64x61  image  from 
4006  to  8192.  While  the  number  of  unknowns  remains  the  same  (4096  complex 
values  or  8192  real  values)  there  are  now  two  projections  that  are  needed  for 
each  row.  Fortunately  the  number  of  operations  remains  the  same. 

The  implementation  of  this  algorithm  was  tested  by  comparing  the  results 
of  the  exact  solution  for  a  small  cylinder.  Using  the  incident  field  as  the  first 
“guess”  for  the  solution  several  iterations  were  calculated  and  the  real  part  of 
the  solution  is  shown  in  Figure  6  26.  This  compares  favorably  with  the  exact 


solution  shown  as  a  dashed  line.  In  this  example  and  the  work  to  follow  one 
iteration  is  defined  as  projecting  a  vector,  x°,  onto  each  of  the  N2  hyperplanes. 

The  speed  at  which  this  algorithm  converges  to  the  correct  answer  is 
shown  in  Figure  6.27.  Here  the  real  part  of  the  field  for  a  section  of  space  near 
the  origin  is  shown  before  iterating  (the  incident  field  is  the  initial  guess)  and 
then  after  one,  two  and  three  iterations.  This  result  is  especially  encouraging 
because  the  first  iteration  has  changed  so  rapidly  towards  the  correct  answer. 

Ramakrishnam  (Ram79j  proposed  that  faster  convergence  of  a  projection 
algorithm  could  be  obtained  by  using  pair  wise  orthogonalization.  As  already 
described,  the  Kaczmarz  approach  will  converge  in  one  step  when  all  the 
equations  are  orthogonal  and  as  the  hyperplanes  become  more  parallel  it  will 
take  longer  for  the  method  to  converge. 

Certainly  the  best  way  to  speed  up  convergence  is  to  first  orthogonalize 
the  system  of  equations.  Then  it  would  be  possible  to  solve  the  system  of 
equations  in  a  single  iteration.  Unfortunately  the  work  required  to 
orthogonalize  the  system  is  identical  to  that  needed  to  find  the  inverse  of  the 
matrix.  In  addition  it  then  would  be  necessary  to  compute  and  store  all  the 
elements  of  the  matrix.  Due  to  the  large  size  of  the  matrix  this  is  not  practical. 

Ramakrishnam  proposed  that  pair  wise  orthogonalization  be  used  to  make 
each  hyperplane  perpendicular  to  the  previous  hyperplane.  With  this  approach 
the  equation  for  the  field  at  each  point  is  made  orthogonal  to  the  preceding  one 
using  the  relation 

A;  =  A-Aj  , 

and 

bj  —  bj  bj_j 

Here  the  new  orthogonalized  system  of  equations  are  denoted  with  the  matrix 
A  and  the  vector  b,  A;  is  one  row  of  the  A  matrix  and  b;  is  the  corresponding 
element  of  the  b  vector. 

While  this  approach  is  not  optimum  it  does  have  the  advantage  that  at 
each  step  storage  for  only  one  extra  row  of  the  A  matrix  is  needed. 
Ramakrishnam  showed  for  a  simple  restoration  problem  that  pair  wise 
orthogonalization  reduced  the  number  of  iteration  needed  to  obtain  a  given 


<Ai,Ai_i> 

<Ai_„Ai_1> 


(6.80) 
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Rea!  Part  of  Scattered  Field 
Radius  =  3X,  Refractive  Index  =  1.01 

Incident  Field  1  Iteration 


Figure  6.27  Three  iterations  of  the  Kaczmarz  algorithm  are  shown  to 
demonstrate  the  convergence  of  the  approach  to  a  single 


mean  squared  error  by  a  factor  of  two.  Since  the  orthogonalization  is  done 
with  several  dot  products  the  reduction  in  iterations  more  than  balances  the 
extra  work. 


An  alternative  approach  is  rearrange  the  order  of  the  equations  to  reduce 
their  interdependency.  This  idea  was  first  published  by  Hounsfield  in  the 
original  patent  for  CT  imaging  [Hou72].  It  seems  reasonable  that  the 
hyperplanes  describing  the  field  at  adjacent  points  (less  than  a  quarter 
wavelength  apart)  would  be  nearly  parallel.  When  the  solution  vector  is  first 
projected  onto  one  plane  then  projecting  onto  a  second  parallel  hyperplane  will 
not  significantly  improve  the  answer.  Thus  convergence  will  be  faster  if  the 
parallel  hyperplane  is  saved  till  later  in  the  iteration  sequence. 

This  idea  is  illustrated  in  Figure  6.28  for  a  two  dimensional  case  with  four 
hyperplanes.  In  each  case  the  order  projections  are  considered  is  indicated  as 
the  hyperplane  number  at  the  end  of  the  line.  First  in  Figure  6.28a  the  two 
sets  of  parallel  planes  are  considered  separately  while  in  Figure  6.28b  the  two 
sets  are  interleaved.  It  is  easy  to  see  that  the  first  ordering  will  take  twice  as 
many  iterations  as  the  second. 

To  calculate  the  scattered  field  using  the  Kaczmarz  algorithm  a  system  of 
equations  is  set  up  that  represents  the  field  at  each  point  as  a  function  of  the 
refractive  index  distribution.  For  ease  of  programming  the  x  vector,  the  field 
at  each  point  in  the  grid,  is  organized  so  that  adjacent  elements  in  the  vector 
represent  the  field  at  adjacent  points  in  the  grid.  Thus  if  the  hyperplanes  are 
considered  in  order  there  will  be  a  high  degree  of  correlation  between  the 
equations  and  convergence  will  be  slow. 

The  degree  of  independence  of  two  equations  can  be  found  by  finding  the 
angle  between  the  two  hyperplanes.  If  two  equations  contain  nearly  the  same 
information  then  the  angle  between  their  respective  hyperplanes  will  be  zero 
while  if  the  two  equations  are  independent  there  will  be  an  angle  of  90  degrees. 
From  standard  vector  theory  the  angle  between  two  hyperplanes  is  defined  as 

COS0  =  ;  J  (6.81) 

\/(<Aj,Aj>  <Aj,Aj>) 

where  A-,  and  Aj  represent  the  rows  of  the  A  matrix  or  the  normal  vector  to  the 
two  hyperplanes. 

The  order  equations  are  considered  is  a  function  of  a  parameter  called 
that  represents  the  change  in  equation  number.  For  each  iteration  of  a  N  x  N 
grid  the  parameter  i  steps  from  0  to  N2-l  and  is  mapped  into  an  equation 
number,  j,  by  the  relation 


The  order  equations  are  considered  can  affect  the  rate  of 
convergence.  If  the  two  sets  of  parallel  lines  are  considered 
separately  (a)  then  the  convergence  is  twice  as  slow  as  it 
would  be  if  they  are  interleaved  (b). 


j  =  (i*^)  mod  N2.  (6.82) 

With  the  proper  choice  of  ^  the  equation  number  j  will  step  through  all  N2 
equations. 

The  table  below  6.2  shows  the  average  and  maximum  cosine  of  the  angle 
for  a  range  of  refractive  indices  between  1.01  and  1.5.  A  system  of  1024 
equations  is  used  to  define  the  scattered  field  over  a  32  x  32  grid.  Thus  when 
the  parameter  ^  is  1  adjacent  equations  are  compared  while  the  equations  are 
as  far  apart  as  possible  when  ^  is  set  to  513. 

Table  6.2.  The  average  and  maximum 
cosine  of  the  angle  between  hyperplanes  is  shown  as  a  function 
of  the  refractive  index  and  the  number  of  equations  skipped. 


Refractive 

Index 

Average 

COS0 

Maximum 

cost? 

1.01 

1 

.000592 

.010902 

513 

.000125 

.001338 

1.1 

1 

.019414 

.1682 

513 

.001981 

.013842 

1.2 

1 

.064287 

.377699 

513 

.005448 

.025288 

1.5 

1 

.238268 

.623253 

513 

.022828 

.06990 

From  the  information  in  this  table  two  points  are  apparent.  First,  as  the 
refractive  index  is  increased  both  the  average  and  the  maximum  cosine  of  the 
angle  increases.  For  small  refractive  indices  the  A  matrix  is  dominated  by  the 
diagonal  terms  and  thus  each  equation  is  nearly  independent  of  the  others. 
This  also  explains  why  the  Kaczmarz  algorithm  converges  much  faster  for 
small  refractive  indices. 

Secondly,  equations  describing  the  field  at  widely  spaced  points  have  a 
large  angle  between  their  respective  hyperplanes.  As  can  be  seen  from  the 
above  table  by  comparing  equations  that  are  widely  separated  the  average  and 
the  maximum  cosine  of  the  angle  are  reduced  by  a  factor  of  10.  Thus  even 


with  a  refractive  index  as  great  as  1.5  the  minimum  angle  between  hyperplanes 
is  increased  to  86  degrees  by  skipping  a  large  number  of  equations. 

The  advantages  gained  by  considering  non  adjacent  equations  is  confirmed 
by  considering  the  convergence  of  the  algorithm.  For  the  system  of  equations 
defined  by  Ax  =  b  the  residue  is  given  by 

Residue  =  <A;,x>-b.  (6.83) 

The  residue  measures  the  distance  between  the  solution  vector  x  and  the 
hyperplane  described  by  A;  and  the  total  residue  is  defined  as  the  sum  of  the 
squares  of  the  residues  from  each  row  of  A.  This  figure  can  then  be  used  as  a 
measure  of  quality  of  the  solution. 

The  total  residue  when  calculating  the  scattered  field  from  a  cylinder  of 
radius  IX  and  refractive  index  1.1  is  shown  in  Figure  6.29.  The  study  was  done 
for  16  iterations  and  compares  the  total  residue  when  the  equations  are 
considered  in  order  is  1)  and  when  the  equations  are  widely  separated  is 
513).  To  reach  any  given  total  residue,  iterating  the  solution  by  considering 
adjacent  equations  takes  twice  as  long  to  converge  as  when  the  equations  are 
widely  separated. 

Considering  widely  separated  points  gives  the  same  benefits  as  pair  wise 
orthogonalization  (convergence  is  twice  as  fast)  but  without  the  extra  work.  In 
addition  when  the  equations  are  widely  separated  the  planes  are  nearly 
perpendicular  and  thus  there  is  little  to  be  gained  by  pair  wise 
orthogonalization. 

It  is  also  possible  to  study  the  effect  of  projecting  the  field  in  a  non 
sequential  fashion  by  considering  the  field  after  one  iteration.  Figure  6.30 
shows  the  exact  (dashed  line)  and  the  Kaczmarz  field  (solid  line)  from  a 
cylinder  of  radius  IX.  The  cylinder  has  a  refractive  index  of  1.01  and  the 
Kaczmarz  iteration  is  carried  out  over  a  32  x  32  grid.  In  addition,  to 
emphasize  the  difference  made  by  changes  in  6%,  the  initial  guess,  x0,  for  the 
field  was  the  incident  field. 

While  at  first  glance  all  of  the  Kaczmar(z  fields  are  very  poor 
approximations  to  the  exact  field  it  is  important  to  note  that  after  one  more 
iteration  all  of  them  have  converged  to  the  correct  answer.  Thus  the  only 
effect  of  altering  the  order  of  the  equations  is  to  change  the  rate  of 
convergence. 

The  reason  for  the  difference  in  one  iteration  of  the  Kaczmarz  is  best  seen 
by  comparing  the  field  when  ^  is  one  and  1023.  Recall  that  skipping  1023 
equations  in  a  system  of  1024  equations  projects  adjacent  equations,  like  the 


^  \  ^ 


original  case  with  ^  =  1,  but  instead  the  equations  are  treated  in  reverse  order. 

Since  the  diagonal  elements  of  the  A  matrix  dominate  the  matrix  with  a 
refractive  index  of  1.01  each  projection  only  modifies  the  field  at  one  point  in 
the  grid.  Thus  with  =  1  the  field  is  calculated  in  the  same  direction  as  the 
incident  field  travels,  from  negative  to  positive  y.  This  means  that  the  field 
inside  the  cylinder  has  been  calculated  before  the  field  is  calculated  at  the 
receiver  line.  On  the  other  hand  when  the  equations  are  projected  with 
=  1023  the  field  at  the  receiver  line  is  calculated  before  the  field  has  been 
adjusted  for  the  scattering  inside  the  cylinder. 

Unlike  the  fixed  point  algorithms  (Bom  and  Rytov  series)  the  Kaczmarz 
algorithm  should  always  converge  to  the  correct  scattered  field.  The  speed  the 
series  converges  is  related  to  the  average  orthogonality  of  the  defining  equations 
and  as  shown  in  Table  6.2  the  equations  become  more  dependent  as  the 
object’s  refractive  index  increases.  Thus  the  Kaczmarz  algorithm  converges 
more  slowly  for  objects  with  large  refractive  index. 

Two  other  factors  determine  the  rate  of  convergence  of  the  Kaczmarz 
series.  As  shown  in  Figure  6.29  the  order  the  equations  are  considered  changes 
the  rate  of  convergence  of  the  series.  For  this  reason  most  of  the  simulations 
shown  in  this  work  were  done  using  a  ^  equal  to  N2/2  +  l  where  the  size  of  the 
grid  is  NxN.  While  it  hasn’t  been  investigated,  even  better  performance  could 
possibly  be  obtained  by  alternating  different  values  for 

Finally  quantization  and  sampling  errors  slow  the  rate  of  convergence. 
This  is  shown  in  Figure  6.31  for  a  sampling  interval  of  .IX  and  .25X.  In  both 
cases  the  radius  of  the  cylinder  is  equal  to  eight  times  the  sampling  interval 
and  the  calculations  were  performed  over  a  32x32  grid.  The  plots  show  the 
real  (solid  line)  and  imaginary  (dashed  line)  components  of  the  scattered  field  as 
calculated  by  the  Kaczmarz  algorithm  and  an  exact  algorithm  based  on  the 
Bessel  function  expansions  described  in  [Wee64  and  Mor68). 

Each  of  the  plots  in  Figure  6.31  shows  the  scattered  field  after  32 
iterations.  While  the  exact  limits  of  the  Kaczmarz  algorithm  are  difficult  to 
define  it  is  possible  to  say  that  the  Kaczmarz  algorithm  has  converged  for 
refractive  indices  up  to  1.4  with  a  sampling  interval  of  .IX  and  up  to  1.2  with  a 
sampling  interval  of  .25X.  Calculating  the  scattered  field  from  objects  with 
larger  refractive  indices  might  be  possible  but  will  need  more  accurate 
implementations. 

All  the  results  shown  here  were  sampled  using  rectangular  basis  functions. 
While  according  to  the  sampling  theorem  a  rectangular  set  of  basis  functions 
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Figure  6.31 


The  exact  scattered  field  and  the  result  of  32  iterations  of 
the  Kaczmarz  algorithm  are  compared  here.  In  each  case  the 
real  component  of  the  field  is  shown  as  a  solid  line  while  the 
imaginary  component  is  shown  as  a  dashed  line. 
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can  represent  the  original  continuous  function  without  errors  more  accurate 
integrations  might  be  possible  using  sine  functions  [Joh83,  Ste81].  This  idea 
has  not  been  explored  here. 
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CHAPTER  7 

HIGHER  ORDER  RECONSTRUCTION  ALGORITHMS 


7.1  Introduction 

Reconstructions  based  on  the  theory  in  the  first  five  chapters  of  this  work 
are  based  on  first  order  approximations  to  the  scattered  field.  In  other  words 
for  both  the  Born  and  the  Rytov  approximations  it  is  necessary  to  assume  that 
the  field  inside  the  object  is  equal  to  the  incident  field  and  then  it  is  possible  to 
derive  a  simple  (linear)  expression  for  the  scattered  field  as  a  function  of  the 
object.  Finding  a  better  estimate  for  the  field  inside  the  object  is  the  central 
problem  in  improving  diffraction  reconstructions. 

The  reconstruction  problem  is  more  difficult  than  the  forward  problem 
discussed  in  Chapter  6  because  now  both  the  object  and  the  field  inside  the 
object  are  unknown.  This  means  that  is  is  necessary  to  design  a  procedure  that 
simultaneously  estimates  both  the  object  and  the  field  inside  the  object.  This 
procedure  is  made  more  difficult  because  the  reconstruction  is  formed  by 
illuminating  the  object  by  a  number  of  different  fields  and  the  exact  field  inside 
the  object  must  be  calculated  for  each  view. 

Three  approaches  to  the  inverse  problem  will  be  described  here.  Most 
general  and  therefore  computationally  most  expensive  is  to  write  a  system  of 
equations  that  describes  both  the  field  and  the  object  and  then  find  the 
solution  vector  that  gives  the  smallest  error.  Unfortunately  the  system  of 
equations  is  non-linear  and  some  sort  of  search  procedure  must  be  used  to  find 
the  best  solution.  This  approach  to  the  problem  was  first  discussed  by  Johnson 
et  al[Joh83j. 

Computationally  less  demanding  solutions  to  the  problem  are  based  on 
iterative  algorithms.  A  perturbational  approach  much  like  that  used  to 
generate  the  Born  or  the  Rytov  series  was  first  proposed  by  Jost  and  Kahn 
[Jos52]  and  later  extended  by  Moses  [Mos5t>]  and  Prosser  [Pro69,  Pro76).  This 
approach  was  first  developed  for  quantum  scattering  problems  but  is  equally 
valid  for  the  electromagnetic  and  acoustic  wave  equations. 


A  hybrid  approach  to  the  problem  was  6rst  proposed  by  Johnson  [Joh83] 
and  is  a  two  step  iteration  procedure.  It  is  based  on  the  idea  that  is  is  only 
necessary  to  calculate  the  object  since  for  any  given  object  and  incident  field  it 
is  possible  to  calculate  the  exact  field  inside  the  object.  Thus  the 
reconstruction  procedure  first  estimates  the  object’s  refractive  index 
distribution  and  then  an  estimate  of  the  field  inside  the  object  can  be 
calculated  using  any  one  of  the  procedures  described  in  Chapter  6.  The  key  to 
this  procedure  is  then  to  calculate  a  better  estimate  of  the  object  given  a  better 
estimate  of  the  actual  field  inside  the  object.  This  approach  will  be  described 
in  section  7.4  as  an  example  of  a  fixed  point  iteration. 

7.2  Non  Linear  Approach 

The  most  general  approach  to  estimate  the  object  given  the  scattered  field 
is  to  define  a  solution  space  that  includes  both  the  refractive  index  of  the 
object  and  the  exact  field  inside  the  object  for  each  of  the  views.  Both  the 
non-linear  equations  and  the  number  of  unknowns  combine  to  make  this  a 
difficult  problem.  The  non-linear  nature  of  the  problem  means  that  a  search 
procedure  must  be  used  to  find  the  best  solution  and  unlike  fixed  point  or 
perturbation  methods  it  is  not  possible  to  say  a  priori  how  fast  the  search  will 
reduce  the  error  or  whether  it  will  ever  converge. 

The  unknowns  in  this  problem  are  defined  over  an  NxN  grid  and  consist 
of  the  object  and  the  exact  scattered  field  from  views.  If  the  number  of 
views,  N^,  is  on  the  same  order  as  N  then  there  are  a  total  of  N3  unknowns  and 
thus  at  least  N3  defining  equations  are  needed  for  a  well  behaved  solution 
[Sar81|. 

For  each  view  measurements  of  the  scattered  field  are  taken  and  this 
defines  N^NM  equations  of  the  form 

=  /ut,^')o(r')gtj"-7')df'  (7.1) 

where  uM  ^  is  the  measured  scattered  field  for  a  view  at  angle  <p  and  position  T. 
Both  u^ff),  the  total  field  inside  the  object  for  the  view  at  angle  4>,  and  off), 
the  object’s  refractive  index,  are  unknown.  Since  the  unknowns  are  calculated 
over  a  discrete  NxN  grid  the  integral  above  becomes  a  summation  or 

i)ofr'  jlgfP-T'  i).  (7.2) 

i 

Unfortunately  the  measured  points  only  contribute  N^,NM  equations 
therefore  there  are  more  unknowns  than  equations  by  a  factor  of  approximately 
N  (again  assuming  that  N^  and  NM  are  on  the  same  order  as  N).  The 


additional  equations  are  defined  by  noting  that  the  field  inside  the  object  must 
also  satisfy  the  Helmholtz  equation.  There  are  N^N2  equations  of  the  form 

ut.^OO  =  +  /ut,^ '  )o(T'  )g(P-7'  )dT'  (7.3) 

or  in  discrete  form 

«t,^(7i)  =  uM(Ti)+T25]ut,^,j)o(f'j)g(rj-r'j).  (7.4) 

Thus  the  combination  of  equations  (7.2)  and  (7.4)  define  a  total  of 
N*Nm  +  N,N2  e(luations  which  must  be  solved  for  both  ut^(F)  and  o(7). 

The  difficulty  caused  by  the  large  number  of  unknowns  is  compounded  by 
the  non  linearity  in  the  equations.  In  both  equations  (7.2)  and  (7.4)  the 
product  of  the  two  unknowns,  the  field  and  the  object,  is  convolved  with  the 
Green’s  function  and  this  product  means  that  the  Kaczmarz  algorithm  as  used 
in  Chapter  6  is  no  longer  applicable. 

The  usual  approach  to  solve  a  system  of  non-linear  equations  is  to  define 
the  error  as  a  function  of  the  difference  between  the  left  and  right  side  of  the 
equations.  An  optimum  solution  is  then  formed  by  a  search  procedure  that 
looks  for  a  minimum  in  the  error  function.  One  implementation  of  this 
algorithm  reported  by  Tracy  et  al  (Tra83)  took  7  hours  of  computer  time  on  a 
small  minicomputer  to  find  the  object  over  an  11x11  grid.  Calculating  the 
object  over  a  larger  grid  (at  least  128x128  is  probably  needed  for  medical 
imaging)  would  be  prohibitively  expensive. 

7.3  Perturbation  Algorithms 

As  already  described  in  Chapter  6  perturbation  algorithms  are  an 
important  technique  for  solving  the  scattering  problem.  This  technique  was 
first  used  to  solve  the  inverse  scattering  equation  by  Jost  and  Kahn  and  the 
general  techniques  are  described  in  more  detail  in  the  books  by  Nayfeh  [Nay73, 
Nay81].  A  discrete  version  of  the  work  by  Jost  and  Kahn  was  first  reported  by 
Devaney  in  [Dcv82j. 

A  perturbational  expansion  of  the  forward  scattering  problem  is  found  by 
letting  the  object  function  be  written  in  terms  of  a  small  perturbation 
parameter  c  or 

o(r)  =  f\(r)  (7.5) 

and  the  total  field  written  as  a  polynomial  in  terms  of  the  same  perturbation 
parameter  or 


Now  both  the  object  and  the  field  are  expressed  as  a  function  of  the  free 
variable  c.  At  first  glance  the  problem  is  made  more  difficult  by  the  addition  of 
the  extra  perturbation  parameter  but  by  gathering  together  the  powers  of  e  the 
problem  can  be  easily  solved. 

The  equation  for  the  total  field 

u(r1  =  UincOf)  +  /u(r '  )o(T'  )g(r-r '  )dr '  (7.7) 

is  now  written  as  a  function  of  c  or 

00  00 

£  £‘Ui(7)  =  uine  +  /  £  )eX(7'  Igf77'  W  ■  (7-8) 

i  =  0  i  =  0 

Rearranging  this  equation  as  a  polynomial  function  of  t  both  the  scattered  field 
and  the  object  satisfy  an  equation  of  the  form 

0  =  [u0(r)-Ui  J  +  (7.9) 

<  [u  if7)—/ u0(r'  )*(?'  )g(T-r'  )dr'  j  + 

^[u.^-ZugC7'  )xcr'  )g(r-T'  )d7']  + 

e3(u1ff)-/UoCF')X(f')gfFT')df']  +  •  •  • 

In  order  for  this  equation  to  be  valid  each  coefficient  of  t  in  the  series 
expansion  must  be  identically  equal  to  zero.  This  requirement  is  all  that  is 

necessary  to  solve  the  more  general  problem  for  u;  as  a  function  of  the 

perturbation  parameter  f,  but  in  the  scattering  problem  only  the  solution  for 
f  =  1  is  interesting.  Therefore  by  setting  the  coefficient  of  each  power  of  c  equal 
to  zero  and  then  setting  the  value  of  t  equal  to  one  the  following  equations 
result 

uo  =  uinc  (71°) 

u,ff)  =  /u0(f'  )X(f'  )gfF-T#  )dr'  (7.11) 

and  in  general 

Ui(r)  =  fuhi(r')x(?'  )g(r-T')drr  i>i.  (7.12) 

This  is  the  same  system  of  equations  defined  as  the  Born  series  in  Chapter  6 
therefore  the  same  conditions  define  the  region  of  convergence  for  this 
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where 


us(?)  -  S(u0)  +  S2(u0)  +  S3(u0)  + 
S'(u0)  =  s(s‘  ’(uq)]. 


(7.18) 


(7.19) 


Using  this  notation  (7.16)  can  be  written 

eip  -  E<‘si(uo)  +  (7  20) 

i 

EE<ifjsj(si(uo))  + 

>  j 

j]y;^i<i<ksl(sj(si(u0)))  +  •  •  • . 

i  i  k 

This  expression  can  be  simplified  even  further  by  denoting  the  iterated  kernel 
Sj(Si(u0))  by  the  expression  Sj;(u0).  Now  the  Born  series  is  written 

=  EVSiK)  +  (7.21) 

i 

EEfi<jsji(uo))  + 

'  j 

EEE^SkjiK)))  +  •  •  •  ■ 
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Notice  that  the  first  summation  above  represents  first  order  scattering  from  a 
number  of  different  objects  while  the  second  set  of  summations  represents  all 
possible  second  order  scatterings  from  the  same  ensemble  of  objects. 

A  series  solution  for  the  object  function  is  found  by  gathering  together  the 
coefficients  of  like  powers  of  t.  This  gives  the  following  polynomial  in  t 

cV'(7)  =  <S,(u0)  +  f2(S2(u0)  +  S1I(u0)]  +  (7.22) 

c3[^3( uo)  S2i(u0)  +S12(u0)  +S1j|(u0)J  +  •  •  •  . 
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Just  as  was  done  in  deriving  the  Born  series  each  power 

of  f  is 
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independent  and  therefore  the  following  equalities  can  be  written 

*/>  \A 

V>  =  S,(u0) 

(7.23) 

$2(uo)  =  _SI|(U0) 

(7.24) 
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®3( uo)  ~  ^2i(u0)  S12(u0)  SU|(u0)  ("-25) 

and  in  general 

Si(«o)=-  E  Siii2...in(u0)  n>2.  (7.26) 

ii  +  i2+  •  •  ■  +in  =  i 

Consider  first  the  equation  for  rjj.  The  solution  for  the  object  functions  is 
only  interesting  when  e  =  l  therefore  V’  =  us  and  the  first  equality  above  can  be 
written 

us(T)  =  / u0(T"  )o,(T'  )g(F-T'  )dr ' .  (7.27) 

This  equation  represents  the  measured  scattered  field  as  the  first  order 
scattered  field  from  the  object  Oj  and  by  the  Fourier  Diffraction  Theorem  this 
equation  can  be  solved  exactly  for  o,.  Because  the  object  is  illuminated  with 
the  incident  field  this  is  true  for  all  experiments,  regardless  of  the  size  of  the 
object  and  its  refractive  index.  This  contrasts  with  first  order  diffraction 
tomography  where  the  object  is  modulated  by  the  total  field  and  thus  the 
Fourier  Diffraction  Tomography  is  only  valid  when  the  total  field  can  be 
approximated  by  the  incident  field. 

While  equation  (7.23)  only  expresses  the  scattered  field  from  one  view  of 
the  object  it  is  possible  to  combine  the  scattered  field  from  a  number  of 
different  views  and  then  use  the  first  order  reconstruction  algorithm  described 
in  Chapter  4.  Thus  the  result  of  the  first  order  reconstruction  algorithms 
described  in  Chapter  4  and  5  is  exactly  equal  to  Oj. 

The  second  order  object  is  slightly  more  difficult  to  compute.  From 
equation  (7.23)  the  first  order  scattering  from  the  second  order  object  is  given 
by 

^2( uo)  ~ -S„(uo).  (7.28) 

The  expression  Su(u0)  represents  the  second  order  scattering  from  the  first 
order  object  and  is  easily  calculated  because  ot  has  already  been  computed. 
The  higher  order  terms  follow  in  a  similar  fashion  except  a  number  of  partial 
scattered  fields  are  summed  to  find  the  first  order  scattered  field  from  O;. 

The  procedure  used  to  calculate  the  higher  order  object  is  slightly  different 
from  that  of  the  first  order  object  because  of  the  location  of  the  receiver  line. 
The  first  order  object  is  a  function  of  the  scattered  field  and  the  placement  of 
the  receiver  line  is  limited  by  experimental  constraints.  On  the  other  hand  the 
higher  order  fields  are  defined  on  a  rectangular  grid  since  an  FFT  based 
implementation  of  the  Born  integral  is  the  most  efficient  procedure.  Thus  for 


the  higher  order  terms  the  field  are  calculated  over  the  entire  grid  and  th'*n 
only  the  field  along  one  side  of  the  grid  is  used  as  input  to  the  reconstruction 
procedure. 

In  summary  the  algorithm  for  reconstructing  the  object  using  the  higher 
order  Born  series  is 

•  Use  the  measured  fields  and  the  first  order  Born  reconstruction 
algorithm  to  compute  Oj. 

•  For  each  i>l  do  the  following. 

•  Calculate  the  higher  order  scattered  field  from  each  of  the 
already  computed  objects  (See  equation  (7.26)). 

•  Use  the  first  order  inversion  algorithm  to  invert  S;  and  find  O;. 

•  Sum  up  each  of  the  O;  to  get  the  object  reconstruction. 

Notice  that  this  algorithm  is  “exact.”  Except  for  the  numerical  approximations 
needed  for  the  reconstruction  procedure  there  are  no  mathematical 
approximations  to  limit  the  quality  of  the  reconstruction. 

The  most  expensive  part  of  this  algorithm  is  not  doing  the  reconstructions 
but  instead  in  computing  the  higher  order  scattered  field,  S;.  While  it  is  easy 
to  implement  a  fast  algorithm  to  compute  each  partial  field  in  S;  the  total 
number  of  integrals  increases  rapidly  with  each  succeeding  iteration.  Table  7.1 
shows  the  number  of  partial  fields  and  the  total  number  of  integrals  needed  for 
each  of  the  first  twenty  iterations.  Since  each  integral  takes  a  constant  amount 
of  CPU  time,  no  matter  how  it  is  implemented,  the  practical  limit  of  this 
algorithm  with  today’s  computers  is  certainly  under  ten  iterations. 

The  convergence  of  this  series  is  dependent  on  the  convergence  of  both  the 
forward  Born  series  shown  in  equation  (7.13)  and  the  object  series  shown  in 
equation  (7.15).  Thus  if  either  series  is  divergent  then  this  reconstruction 
procedure  will  also  diverge  and  produce  an  undefined  answer. 

The  convergence  of  the  forward  series  was  discussed  in  Chapter  6.  This, 
for  example,  showed  that  for  an  object  of  radius  2\  the  scattered  fields  could 
be  calculated  using  the  Born  series  only  for  objects  with  a  refractive  index  of 
less  than  about  11%.  This  puts  a  severe  limitation  on  the  type  of  objects  that 
can  be  reconstructed  with  the  higher  order  Born  series.  For  objects  that  do  fall 
within  the  allowable  range  the  reconstruction  with  the  higher  order  Born  series 
should  be  quantitatively  more  accurate  than  that  done  with  a  first  order 
algorithm. 


Table  7.1.  The  number  of  partial  field 
terms  and  integrals  needed  to  calculate  each 
iteration  of  the  inverse  higher  order  Born  series. 


Iteration 

Terms 

Integrals 

1 

0 

0 

2 

1 

2 

3 

3 

7 

4 

7 

19 

5 

15 

47 

6 

31 

111 

7 

63 

255 

8 

127 

575 

0 

255 

1279 

10 

511 

2815 

11 

1023 

6143 

12 

2047 

13311 

13 

4095 

28671 

14 

8191 

61439 

15 

16383 

131071 

16 

32767 

278527 

17 

65535 

589823 

18 

131071 

1245183 

19 

262143 

2621439 

•STww 
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The  seriousness  of  this  limitation  is  further  seen  by  recalling  that  the 
higher  order  Born  series  only  converges  when  the  first  order  field  closely 
approximates  the  total  scattered  field.  This  is  the  same  condition  that 
determines  the  accuracy  of  the  first  order  reconstruction  algorithms  so  the 
ultimate  improvement  is  limited  by  the  quality  of  first  order  reconstructions. 
Thus  it  will  not  be  possible  to  image  any  object  with  a  larger  refractive  index 
or  radius  than  those  in  Figure  7.1  using  an  algorithm  based  on  the  Born  series. 

The  convergence  of  the  object  series  has  yet  to  be  determined.  Jost  and 
Kahn  in  [Jos52]  report  the  convergence  of  the  higher  order  Born  inversion 
procedure  for  two  quantum  mechanical  scattering  experiments. 

While  the  above  derivation  of  the  perturbation  approach  has  expanded  the 
object  as  a  perturbation  about  zero  it  is  also  possible  to  consider  the  object  to 
be  a  small  perturbation  of  a  known  object.  This  generalization  is  known  as  the 
Distorted  Wave  Born  Approximation  (DWBA)  and  is  described  in  [Tay83, 
New66,  Dev83  and  Bey85j.  This  procedure  is  made  more  difficult  because  the 
Green’s  function  used  in  the  integral  now  represents  the  scattered  field  from  a 
point  source  with  the  effect  of  the  known  object  included.  The  convergence  of 
this  procedure  is  not  known. 

7.4  Fixed  Point  Algorithms 

A  third  untested  approach  to  solve  the  inverse  diffraction  problem  is  to  use 
a  fixed  point  algorithm.  While  the  overall  algorithm  is  relatively 
straightforward  it  is  necessary  to  perform  a  first  order  reconstruction  of  the 
object  when  illuminated  by  an  arbitrary  field.  This  is  much  more  difficult  than 
the  first  order  reconstruction  algorithms  based  on  plane  wave  illumination 
described  in  Chapter  4.  (The  synthetic  aperture  approach  does  use  point 
sources,  but  since  a  different  phase  is  added  to  the  scattered  field  for  each 
transmitter  position  a  plane  wave  is  synthesized.) 

A  fixed  point  algorithm  for  calculating  the  object  that  scattered  a 
measured  field  is  based  on  the  equation 

o  =  f(o)  (7.29) 

where  o  is  the  desired  object  function.  Within  the  limits  of  convergence  of  the 
series,  an  initial  guess  Oj_j  can  be  improved  upon  by  iterating  the  equation 

oj  =  f(Oi_,).  (7.30) 

The  exact  form  of  the  iteration  function,  f,  can  take  a  number  of  different 
forms.  Given  only  the  scattered  fields  then  either  the  Born  or  the  Rytov 
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approximation  can  be  used  to  make  an  initial  guess  for  the  object.  Both  of 
these  approximations  assume  that  the  field  inside  the  object  can  be 
approximated  by  the  incident  field  and  this  is  the  major  source  of  error  in  the 
first  order  reconstruction  procedures. 

It  seems  reasonable  that  a  better  estimate  of  the  object  function  could  be 
found  if  the  field  inside  the  object  is  known.  While  it  is  not  possible  to  know 
this  without  knowing  the  object  first,  a  good  estimate  of  the  field  should  be 
possible  given  an  estimate  of  the  object.  If  the  initial  estimate  of  the  object  is 
“good  enough”  then  a  calculation  of  the  scattered  field  given  this  estimate  of 
the  object  should  be  more  accurate  then  just  using  the  incident  field. 

Using  this  new  estimate  for  the  field  inside  the  object  a  better  estimate  of 
the  object  function  should  be  possible.  The  general  iteration  formula  can  now 
be  written 

Oj  =  Reconstruct(Estimate  Fieldfoj.j))  (7.31 ) 

Here  the  function  labeled  “Estimate  Field”  consists  of  estimating  the  total 
fields  inside  the  object  o^,  and  the  function  labeled  "Reconstruct”  consists  of 
estimating  the  object  given  the  total  field  in  the  object.  The  first  iteration  is 
the  simplest  since  the  estimate  of  the  total  field  in  the  object  is  simply  the 
incident  field.  Higher  order  iterations  are  made  even  more  difficult  since  it  is 
necessary  to  estimate  the  fields  inside  the  object  for  each  of  the  views. 

Any  number  of  means  can  be  used  to  implement  the  two  different  steps  in 
the  algorithm.  Estimating  the  field  inside  the  object  can  be  done  with  any  of 
the  procedures  described  in  Chapter  6,  including  the  Bom  and  Rytov  series  or 
the  algebraic  approach.  The  procedure  to  use  would  depend  on  whether  the 
object  falls  in  the  algorithm’s  region  of  convergence  and  on  the  efficiency  of  the 
algorithm  with  the  available  hardware. 

Inverting  the  total  fields  to  get  an  estimate  for  the  object  is  the  most 
difficult  part  of  the  algorithm.  The  Fourier  Diffraction  Theorem  only  applies  to 
objects  illuminated  with  a  plane  wave  so  a  more  general  approach  is  needed. 
One  solution  to  this  problem  was  proposed  by  Vezzetti  and  Aks  [Vez79],  Their 
work  still  assumes  plane  waves  inside  the  object  but  now  the  field  inside  the 


object  is  modified  by  the  average  refractive  index.  With  this  approach  they  do 
show  an  improvement  in  the  quality  of  the  reconstruction  but  it  is  doubtful 
whether  this  approach  would  be  accurate  enough  for  a  fixed  point  algorithm. 

A  complete  solution  for  the  object  given  an  arbitrary  set  of  illuminating 
fields  would  undoubtably  be  based  on  a  least  squared  approach.  While  there 
are  enough  equations,  given  the  field  everywhere,  to  determine  the  object  the 


system  of  equations  would  be  very  unstable  because  the  Green’s  function  only 
samples  one  arc  of  the  scattering  potential’s  Fourier  transform.  This  means 
that  if  the  field  is  known  over  an  NxN  grid  and  it  is  desired  to  calculate  the 
object  over  the  same  grid  then  there  would  be  a  total  of  N2N^  equations 
defining  the  N2  unknowns.  Thus  the  system  of  equations  to  determine  the 
object  is  overdetermined  and  any  error  in  the  field  estimates  will  lead  to  an 
inconsistent  set  of  equations.  A  least  squared  approach  could  then  be  used  to 
find  the  solution  vector  that  best  satisfies  the  defining  equations. 

The  convergence  of  this  method  is  unknown.  Like  the  fixed  point 
methods  discussed  in  Chapter  6  it  is  necessary  for  the  “derivative”  of  the 
function  be  less  than  one  in  some  region  for  the  algorithm  to  converge.  It 
probably  isn’t  unreasonable  to  assume  that  this  method  will  only  converge 
when  the  first  order  estimate  of  the  object  is  “good.”. 
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CHAPTER  8 
CONCLUSIONS 


A  number  of  ideas  are  new  to  this  work.  While  they  have  served  to  answer 
a  number  of  questions  about  diffraction  tomography  there  remains  much  work 
to  be  done.  This  chapter,  therefore,  reviews  the  state  of  the  art  of  diffraction 
tomography  as  presented  by  this  work  and  indicates  directions  for  future 
research. 

Chapter  2  reviewed  the  wave  equation  and  its  integral  solution.  While 
this  material  is  well  known  among  people  doing  research  in  diffraction 
tomography  and  inverse  scattering  its  presentation  here  emphasized  the 
common  mathematical  problems  in  acoustic  and  electromagnetic  scattering. 
For  this  reason  all  distances  were  expressed  in  wavelengths  and  the  object 
function  represented  the  (complex)  refractive  index  variation  of  an 
inhomogeneity  for  either  acoustic  or  electromagnetic  waves.  Researchers  more 
concerned  with  experimental  work  will  want  to  use  the  relationships  presented 
in  Chapter  2  to  convert  the  results  presented  in  the  remainder  of  this  work  to 
more  physical  quantities. 

Finally  Chapter  2  also  presented  two  approximations,  the  Born  and  the 
Rytov,  which  allow  linearized  versions  of  the  wave  equation  to  be  written. 
These  two  first  order  perturbational  approximations  are  important  because 
they  allow  simple  inversion  algorithms  to  be  derived.  Since  these 
approximations  are  so  critical  to  first  order  diffraction  tomography  the 
mathematical  limitations  of  each  approximation  are  also  discussed. 

The  Fourier  Diffraction  Theorem  relates  the  scattered  field  measured  on  a 
line  to  the  Fourier  transform  of  the  object  and  is  presented  in  Chapter  3.  This 
theorem  is  only  true  when  either  the  Born  or  the  Rytov  approximation  is  valid 
but  it  has  generated  much  excitement  in  the  research  community. 

The  Fourier  Diffraction  Theorem  was  derived  by  two  different  methods  in 
this  work.  Both  approaches  to  the  Fourier  Diffraction  Theorem  lead  to  the 
same  relationship  between  the  scattered  field  and  the  object’s  Fourier 
transform.  The  conventional  approach  is  to  decompose  the  Green’s  function, 
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the  field  scattered  by  a  point  scatterer,  into  plane  waves  and  simply  substitute 
this  result  into  the  integral  solution  to  the  wave  equation.  A  second,  new 
approach,  is  to  consider  the  Fourier  Diffraction  Theorem  entirely  in  the  Fourier 
domain.  This  method  points  toward  a  more  natural  computer  implementation 
and  was  exploited  in  Chapter  6  for  computing  better  approximations  to  the 
scattered  field. 

The  remainder  of  Chapter  3  discussed  experimental  methods  for  collecting 
enough  scattered  data  so  that  a  unique  estimate  for  the  object  can  be  formed. 
The  potential  methods  described  include  using  a  plane  wave  to  illuminate  the 
object,  synthesizing  plane  waves  much  like  what  is  done  in  phase  array  antenna 
design  and  broadband  (in  time)  incident  fields. 

Chapter  4  discussed  two  mathematical  algorithms  for  inverting  the 
scattered  data  to  estimate  the  object’s  (complex)  refractive  index.  Much  like 
conventional  (straight  ray)  tomography  there  are  two  approaches  that  can  be 
used  to  invert  the  scattered  data.  These  two  algorithms,  often  described  as 
interpolation  in  the  space  domain  and  frequency  domain,  were  presented  here 
and  their  algorithmic  complexity  was  discussed. 

In  addition  several  signal  processing  concerns  were  examined  in  Chapter  4. 
By  calculating  the  Mean  Squared  Error  between  the  object  and  the 
reconstruction  it  was  concluded  that  zero  padding  each  projection  is  a  good 
way  to  reduce  the  interpolation  error  in  the  frequency  domain.  On  the  other 
hand,  using  a  Hamming  window  to  shape  the  projection  and  reduce  the  effect 
of  the  finite  aperture  severely  attenuates  the  high  frequency  information  in  the 
projection  and  increases  the  error. 

The  limitations  of  first  order  diffraction  tomography  were  discussed  in 
Chapter  5.  Two  types  of  errors  limit  the  quality  of  the  reconstruction: 
mathematical  limitations  caused  by  the  approximations  used  to  derive  the 
Fourier  Diffraction  Theorem  and  experimental  limitations  caused  by  the  ability 
to  only  collect  a  finite  amount  of  data.  The  mathematical  limitations  are  the 
most  severe.  In  deriving  the  Born  and  the  Rytov  approximations  it  was 
necessary  to  assume  that  the  scattered  fields  were  small  compared  to  incident 
fields.  This  is  equivalent  to  saying  that  the  object  must  be  weakly  scattering 
for  the  first  order  diffraction  tomography  algorithms  to  hold  and  if  this 
condition  isn’t  met  then  the  reconstruction  will  have  serious  artifacts. 

The  limits  of  first  order  diffraction  tomography  are  easily  described  in 
terms  of  the  magnitudes  of  the  scattered  fields  but  a  more  meaningful  measure 
is  to  study  the  range  of  objects  where  the  approximations  are  valid.  This  was 
done  in  Chapter  5  by  calculating  the  exact  scattered  fields  from  a  large  number 
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of  cylinders  and  then  making  an  estimate  of  the  object  assuming  that  the  first 
order  diffraction  tomography  algorithms  are  valid.  Thus  it  was  concluded  that 
the  Born  approximation  is  valid  when  the  product  of  the  diameter  of  the 
cylinder  (in  wavelengths)  and  the  absolute  value  of  the  refractive  index  change 
is  less  than  0.5.  On  the  other  hand  the  size  of  the  object  is  not  as  critical  to 
the  Rytov  approximation.  Instead  the  refractive  index  change  is  the  limiting 
factor  and  reconstructions  based  on  the  Rytov  approximation  are  good  as  long 
as  the  refractive  index  of  the  object  is  less  than  a  few  percent. 

The  experimental  limitations,  on  the  other  hand,  can  always  be  minimized 
by  collecting  more  data.  Thus  it  is  clear  that  interpolation  error  can  always  be 
reduced  by  increasing  the  number  of  projections  or  the  number  of  samples  per 
projection.  Another,  less  obvious,  limitation  is  the  finite  aperture  of  the 
projection.  Unlike  conventional  (straight  ray)  tomography  where  the  projection 
of  a  finite  sized  object  has  a  finite  length,  the  same  is  not  true  for  scattered 
fields.  With  diffraction  tomography  the  scattered  field  never  goes  to  zero  and 
the  sampling  interval  for  the  projection  must  be  carefully  balanced  to  prevent 
aliasing  but  yet  large  enough  to  measure  the  high  frequency  information  far 
from  the  center  of  the  projection.  An  expression  for  this  relationship  was 
derived  in  Chapter  5  and  several  reconstructions  were  presented  with  different 
sampling  intervals  to  confirm  the  optimum  sampling  interval. 

The  limitations  of  first  order  reconstruction  algorithms  were  addressed  in 
Chapters  6  and  7.  The  most  severe  limitation  is  caused  by  the  first  order 
perturbation  models  assumed  in  deriving  the  Fourier  Diffraction  Theorem. 
Thus  Chapter  6  discussed  several  approaches  to  more  accurately  model  the 
scattered  fields.  With  one  of  these  more  accurate  models  it  should  then  be 
possible  to  invert  the  relationship  and  derive  a  better  reconstruction  algorithm. 
A  survey  of  several  possible  approaches  to  inverting  the  scattered  data  is 
presented  in  Chapter  7. 

Since  better  reconstructions  will  be  based  on  more  accurate  models  of  the 
field  inside  the  object  two  approaches  to  more  accurately  model  the  scattered 
field  were  presented  in  Chapter  6.  The  most  severe  limitation  of  first  order 
algorithms  is  the  assumption  that  the  field  inside  the  object  is  approximately 
equal  to  the  incident  field.  Thus  when  this  condition  is  not  valid  the  Born  and 
the  Rytov  approximations  are  no  longer  valid. 

The  simplest  technique  is  to  assume  that  the  perturbational  model  used  to 
derive  the  Fourier  Diffraction  Theorem  is  approximately  correct  and  simply 
include  more  of  the  higher  order  terms  The  result  is  a  series  of  terms  much 


like  a  Taylor  series.  This  is  an  iterative  procedure  and  was  applied  to  both  the 
Born  and  the  Rytov  approximations. 

An  important  measure  of  any  series  is  a  description  of  its  region  of 
convergence.  In  this  case  the  region  of  convergence  is  a  function  of  the  entire 
object  and  the  results  presented  in  Chapter  6  were  simplified  by  considering  the 
convergence  as  a  function  of  size  and  refractive  index  of  simple  objects.  Thus 
the  region  of  convergence  can  be  described  by  two  parameters  and  all  objects 
outside  this  region  (because  they  are  larger  or  have  a  greater  refractive  index 
change)  will  cause  the  series  to  diverge. 

The  series  described  in  Chapter  6  were  calculated  by  sampling  the  object 
and  the  fields  and  then  using  an  efficient  algorithm  based  on  Fourier 
transforms.  In  each  case  the  scattered  field  was  calculated  by  multiplying  a 
function  of  the  object  by  a  field  and  then  convolving  this  “scattering  potential’’ 
with  the  Green’s  function.  The  convolution  represents  the  most  expensive  part 
of  the  algorithm  and  can  be  efficiently  calculated  using  FFT’s. 

The  convergence  properties  of  the  Born  and  the  Rytov  series  were 
determined  by  a  binary  search  procedure.  Thus  for  a  given  size  the  refractive 
index  of  the  object  was  varied  till  a  point  was  found  were  the  series  converged 
for  all  refractive  indices  that  were  smaller  and  diverged  if  the  refractive  index 
was  larger  than  this  point.  By  varying  the  size  of  the  object  it  was  possible  to 
show  a  region  of  convergence  as  a  function  of  both  object  size  and  refractive 
index. 

The  simulations  of  the  Born  series  showed  it  to  converge  only  when  the 
first  iteration  is  an  accurate  estimate  of  the  field  inside  the  object.  Thus  the 
phase  change  of  the  field  as  it  travels  through  the  object  is  a  good  indication  of 
not  only  the  quality  of  a  first  order  Born  reconstruction  but  also  describes  the 
region  of  convergence  of  the  Born  series. 

The  convergence  of  the  Rytov  series  is  more  surprising.  For  all  cases 
studied  the  Rytov  series'  region  of  convergence  includes  the  region  of 
convergence  for  the  Bom.  This  is  especially  surprising  since  the  first  order 
Born  and  Rytov  approximations  have  different  regions  of  validity. 

In  addition  Chapter  6  also  presented  a  study  of  the  effects  of  attenuation 
on  both  the  Born  and  the  Rytov  series.  A  key  part  of  this  work  is  the  idea 
that  attenuating  plane  waves  can  be  described  either  in  terms  of  a  solution  to 
the  wave  equation  or  in  the  Fourier  domain  In  a  non-attenuating  media  the 
two  approaches  are  identical  since  plane  wave  solutions  to  the  wave  equation 
are  also  Fourier  waves. 


Considering  an  attenuating  plane  wave  in  the  Fourier  domain  makes  it 
possible  to  calculate  the  higher  order  Born  and  Rytov  series  for  attenuating 
media.  While  the  algorithms  remain  the  same  there  is  a  significant  difference 
in  its  convergence  properties.  Since  the  energy  in  the  field  is  attenuated  as  it 
travels  away  from  the  scattering  site  the  region  of  convergence  for  both  the 
Born  and  the  Rytov  series  is  increased  as  the  attenuation  of  the  media  is 
increased.  Thus  the  attenuation  of  the  media  balances  the  extra  field  caused 
by  a  larger  scattering  potential. 

A  second  approach  to  calculating  the  scattered  fields  from  a  known  object 
was  also  discussed  in  Chapter  6.  By  sampling  the  object  and  the  fields  a  set  of 
discrete  equations  can  be  written  that  relate  the  field  and  the  object.  Without 
using  any  approximations  it  is  then  possible  to  express  the  field  as  the  solution 
of  a  linear  matrix  equation. 

While  the  form  of  the  matrix  equation  is  simple,  the  large  amount  of  data 
makes  this  problem  difficult  to  compute  directly  with  today’s  computers 
Instead  it  was  necessary  to  use  an  iterative  technique  known  as  the  Kaczmarz 
approach  to  solve  the  matrix.  While  the  iterative  technique  used  can  be  shown 
theoretically  to  always  converge,  numerical  errors  limit  the  range  of  objects  to 
those  that  have  a  refractive  index  change  of  less  than  2CM0ro. 

The  rate  of  convergence  of  this  method  is  only  a  function  of  the 
orthogonality  of  the  defining  equations.  Thus  when  the  object  has  a  small 
refractive  index  the  defining  equations  are  nearly  orthogonal  and  the  Kaczmarz 
approach  quickly  converges  to  the  correct  field.  On  the  other  hand  as  the 
refractive  index  is  increased  the  hyperplanes  defined  by  the  equations  become 
nearly  parallel  and  convergence  is  much  slower.  Since  the  Kaczmarz  approach 
treats  each  equation  for  the  field  separately  faster  convergence  is  often  possible 
by  sequencing  the  equations  so  that  each  equation  is  nearly  parallel  to  the  one 
before  it. 

Finally  Chapter  7  presented  a  survey  of  several  techniques  for 
reconstructing  an  object  without  using  first  order  approximations.  The  most 
difficult  part  of  this  problem  is  that  it  now  necessary  to  actually  compute  the 
field  inside  the  object.  In  first  order  diffraction  tomography  the  field  inside  the 
object  is  assumed  to  be  a  plane  wave  but  this  can’t  be  true  with  higher  order 
approximations.  Since  it  is  necessary  to  illuminate  the  object  from  a  number 
of  different  directions  to  perform  the  reconstruction  a  calculation  of  the  field  is 
necessary  for  each  view.  The  large  number  of  equations  makes  this  a  difficult 
and  expensive  process. 


A  straightforward  approach  is  to  write  a  system  of  equations  that 
describes  both  the  field  inside  the  object  and  the  refractive  index  of  the  object. 
It  should  then  be  possible  to  solve  this  system  of  equations  for  both  the  field 
and  the  object.  Unfortunately  it  is  a  non-linear  system  of  equations  because 
the  defining  equations  are  a  function  of  the  product  of  the  two  unknowns.  For 
this  reason  it  is  necessary  to  use  some  type  of  search  procedure  to  solve  for 
both  the  fields  and  the  object. 

A  second  approach,  first  used  in  high  energy  physics  and  described  in 
Chapter  7,  is  to  do  a  perturbation  expansion  for  the  object.  This  is  similar  to 
the  Born  and  the  Rytov  series  described  in  Chapter  6  but  now  the  object  is 
assumed  to  consist  of  a  series  of  components. 
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The  convergence  of  this  approach  is  a  function  of  two  series.  Since  this 
approach  is  based  on  a  Born  series  expansion  for  the  scattered  field  it  is  only 
valid  when  the  field  inside  the  object  can  be  described  as  a  converging  Born 
series.  As  seen  in  Chapter  6  this  is  a  rather  severe  limitation.  In  addition  the 
object  is  expressed  as  a  separate  series  expansion  and  for  this  approach  to 
converge  it  is  necessary  for  both  the  Born  series  and  the  object  series  to 
converge. 

Finally  a  third  approach,  described  in  Chapter  7,  is  to  make  a  first  order 
estimate  for  the  object  and  then  use  this  object  to  calculate  a  better  estimate 
for  the  field  insiue  the  object.  Like  the  Born  and  the  Rytov  series  described  in 
Chapter  6  this  is  a  fixed  point  algorithm.  This  approach  is  made  even  more 
difficult  than  first  order  reconstruction  algorithms  since  it  is  necessary  to 
calculate  an  estimate  of  the  object  given  an  arbitrary  illuminating  field.  Since 
each  projection  is  no  longer  independent  the  Fourier  Diffraction  Theorem  is  not 
valid  and  a  reconstruction  procedure  will  need  to  look  at  all  the  scattered  data 
simultaneously.  This  can  be  easily  done  using  a  matrix  formulation  but  there 
is  a  severe  performance  penalty. 

The  convergence  properties  of  this  particular  series  is  not  known  although 
it  is  probably  reasonable  to  assume  that  the  region  of  convergence  will  be  a 
function  of  the  quality  of  the  first  order  estimate  of  the  field.  If  using  the  first 
order  estimate  of  the  field  is  not  better  than  the  original  assumption  to  use  the 
incident  field  then  certainly  the  series  will  diverge.  This  condition  represents  a 
severe  limitation  for  the  technique. 

Future  work  on  this  problem  could  continue  in  several  areas.  The 
perturbational  approach  has  a  limited  range  of  convergence  but  for  objects  that 
fall  within  this  range  a  better  quantitative  estimate  of  the  object  should  be 
possible  than  that  which  is  possible  using  first  order  algorithms.  The  same  is 


also  true  for  the  fixed  point  approach  but  more  work  is  needed  to  determine 
the  range  of  convergence. 

Certainly  the  only  guaranteed  approach  to  solve  the  inverse  scattering 
problem  is  to  find  a  solution  to  a  non-linear  set  of  equations.  There  are  a 
number  of  algorithms  that  can  be  used  but  the  large  number  of  equations  (a 
128x128  reconstruction  has  over  2  million  unknowns  and  equations)  makes  this 
a  very  difficult  problem. 
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